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SUMMARY 


C 


c 


This  report  summarizes  the  second  year  work  under  the  general  heading 
of  "Effect  of  Local  Material  Imperf ections  on  the  Buckling  Behavior  of 
Composite  Structural  Elements". 

It  describes  two  important  areas:  one  of  delamination  buckling  of 
cylindrical  shells  under  axial  compression,  and  second  buckling  of  mulr.i*- 
annular  plates,  which  emphasizes  the  effect  of  holes  and  foreign  rigid 
Inclusions. 

The  findings  of  the  research  are  reported  in  two  parts.  The  first  part 
(Part  A)  deals  with  buckling  of  axially'-loaded  delaminated  thin  cylindrical 
shells.  The  geometry  is  assumed  to  be  virtually  isotropic  and  the  size  and 
position  of  the  delamination,  which  extends  around  the  entire  circumference 
of  the  cylinder,  are  varied.  Moreover,  two  sets  of  boundary  conditions  are 
considered  simply  supported  and  clamped. 

The  second  part  deals  with  buckling  of  multi.-annular  plates,  made  out 
of  different  materials,  supported  in  various  ways,  and  subjected  to  radial 
compressive  loads.  The  rigidities  of  the  parts  are  varied  in  order  for  the 
analysis  to  include  plates  with  holes,  rigid  inclusions,  and  ring 
stiffeners. 
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INTRODUCTION 

The  use  of  f iber-reinforced  composite  laminates  is  rapidly 
increasing  due  to  their  excellent  undirectional  properties  and  their 
high  strength  to  weight  ratio.  Recently,  composite  laminates  have  been 
used  in  space  vehicles  in  the  form  of  circular  cylindrical  shells,  as  a 
primary  load  carrying  structure.  In  order  to  use  laminated  composite 
structures  properly  and  effectively,  it  is  important  to  predict  and 
understand  their  behavior  under  the  different  external  actions. 

Imperfections  have  a  very  significant  effect  on  the  structural 
response  and  on  its  load  carrying-  capaci ty.  In  addition  to  the  usual 
geometric  imperfections  that  are  inherent  in  metallic  plates  and  shells, 
laminated  geometries  are  prone  to  such  defects  as  broken  fibers,  cracks 
in  the  matrix  material,  delaminated  areas  (  separatioin  of  adjoining 
plies),  as  well  as  holes  and  small  voids. 

Delamination  is  one  of  the  most  common  failure  modes  in  composite 
materials.  Delamination  is  developed  as  a  result  of  imperfections  in 
the  production  technology  or  due  to  the  effect  of  certain  factors  during 
the  operational  life  of  the  laminate,  such  as  impact  by  foreign  objects. 
The  presence  of  delamination  in  a  laminated  composite  material  may  cause 
local  buckling  and  reduction  in  the  overall  stiffness  of  the  structure, 
which  may  lead  to  early  failure. 


( 


( 


The  problems  of  delaminated  composite  structures  have  received 
attention  in  recent  years.  As  far  as  the  buckling  response  of 
delaminated  structures  is  concerned,  an  extensive  review  is  <-:ven  in 


Ref . [1] . 

The  problem  of  buckling  of  delaminated  cylindrical  shells  has  not 
received  the  deserved  attention.  Very  few  investigations  have  been 
carried  out  in  this  area.  Kulkarni  and  Frederick  [2]  use  a  "branched 
integration"  technique  to  solve  the  problem  of  buckling  of  a  two-layered 
cylindrical  shel 1 .partial ly  debonded,  and  subjected  to  axial 
compression.  They  [2]  consider  the  case  where  the  delamination 
originated  at  the  boundary.  Results  are  reported  [2]  for  different 
lengths  of  debonding  and  inner  to  outer  layer  thickness  ratios.  A 
significant  decrease  in  the  critical  load  is  observed.  The  buckling  of 
stiffened  circular  cylindrical  shells,  with  two  unbonded  orthotropic 
layers,  is  reported  by  Jones  [3].  Jones  [3]  assumes  that  the  layers  do 
not  seperate  during  buckling,  i.e.  the  deformation  of  both  layers  are 
assumed  to  be  the  same.  He  also  examines  the  case  when  one  of  the  two 
unbonded  orthotropic  layers  is  circumferentially  cracked.  His  results 
for  a  cylindrical  shell  made  of  aluminum  with  ablative  outer  layer  and 
subjected  to  hydrostatic  pressure  show  that  the  ablative  layer  had  to  be 
increased  by  50^  in  thickness  in  the  damaged  (debonded)  cylindrical 
shell  in  order  to  obtain  the  same  buckling  load  as  that  of  the  perfect 
cylindrical  shell.  Troshin  [4]  examines  the  effect  of  longitudinal 
delamination  in  a  laminar  cylindrical  shell  on  the  critical  external 
pressure.  The  delamination  is  assumed  to  extend  over  the  entire  length 
of  the  shell.  He  [4],  reports  on  the  effect  of  the  length  and  position 
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CHAPTER  I  I 


BUCKLING  OF  DELAMINATED  CYLINDRICAL  SHELLS 


2.1  MATHEMATICAL  FORMULATION 


effect  of  circumferential  delamination  on  the  critical  loads. 


for  a  laminated  circular  cylindrical  shell  and  subjected  to  compressive 
axial  loading,  is  studied.  The  material  behavior  of  the  cylindrical 
shell  is  assumed  to  be  linearly  elastic.  Delamination  is  assumed  to 
exist  before  the  load  is  applied,  and  extends  along  the  entire 
circumference  of  the  cylindrical  shell,  on  a  surface  parallel  to  the 
reference  surface  (see  Figure  1).  The  location  and  s i 2e  of  the 
delamination  is  arbitrary  and  the  boundary  conditions  are  either  simply 
supported  or  clamped.  Del  ami  nation  separates  the  cylindrical  shell  into 
four  regions  (four  thin-walled  cylindrical  shells),  such  that  each 
region  is  symmetric  with  respect  to  its  own  mid-surface. 

The  governing  equations  for  the  buckling  of  the  delaminated 
cylindrical  shell  are  derived  based  on  the  thin  shell  theory 
(Ki rchhof f-Love  hypotheses),  Donel 1-type  of  kinematic  relations,  and  the 
assumption  of  existence  of  a  membrane  state,  before  buckling  (classical 
approach)  . 


The  geometry  and  sign  convention  are  shown  in  figures  1  and  2 

The  coordinate  system  is  such  that  x  is  measured  from  the  -eft  end. 

Moreover  u  ,v  and  w  (i*l,2,3»M  denote  the  ax i a  1 , c i rcumf erent i a  1 
i  i  i 

and  radial  displacements  of  the  material  points  on  the  mid  surface  of 
each  reg i on, respect ive ly  . 


2.1.1  Kinematic  Relations 

The  axial  and  c i rcumf erent i a  1  displacements  of  any  material  point 

(x , 2  )  are  given  by 

i 


^(x.y.z.)  -  u°(x,y)  -  zL  t# 
vi(x>y*2i)  *  v°(*,y)  -  Zj 


(1) 


where  u°  and  v°  are  components  of  "middle  surface"  material  points  , 

2^  is  measured  from  the  mid  surface  of  each  region  and  the  comma  denotes 
partial  differentiation  with  respect  to  the  index  that  follows. 

The  kinematic  (Oonnel 1 -type)  relations  are  given  by  : 


where  the  superscript  "o"  denotes  reference  strains  and  the  k's  denote 
the  reference  surface  changes  in  curvature  and  torsion. 

Note  that  for  each  region  the  reference  surface  has  been  piaced  at  the 
middle  surface  .  The  reference  surface  strains  are  defined  in  terms  of 
the  displacements  as: 


xx. 


u.  +  w. 
i,x  i.x 


I  -  V,  -  w. /R  +  kv. 

yy<  i,y  j-  t.y 


(3) 


*yi 


U  +  V.  +  W.  W. 

i.y  t,x  L,X  t.y 


and  the  changes  in  curvature  and  torsion  are  given  in  terms  of  the 
transverse  displacement  w  ,by: 


2.1.2  $tress-Strai n  Relations 

Each  lamina  of  Che  composite  shell  can  be  considered  as  orthotropic 
with  two  principal  material  directions  (1  and  2)  parallel  and 
perpendicular  to  the  direct'on  of  the  filaments. 

The  stress-strai n  relations  for  each  lamina  are: 


The  relation  between  the  elements  of  the  stiffness  matrix  [Qj  and  the 
engineering  material  constants  (E^i  .^22  •  c]2  ,v12^  *are  : 


Q11  *  Hi'  (1'V12V21)  •*  Q12  m  V21S11  ' 

(5) 

Q22  “  E227  (1-'J12V21)  5  *33  *  G12 


For  a  general  lamina  (k)  with  fibers  making  an  angle  9^  with  the  x-axis 
then. 
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where 


ft]  *  CT]"1  [Q]  [T]  (7) 

where  [T]  is  the  transformati on  matrix  [5]  . 

It  is  more  convenient  to  deal  with  stress  and  moment  resultants 
rather  than  dealing  with  stresses  in  deriving  the  governing  equations. 
The  stress  and  moment  resultants  are  defined  by 


1Z 
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i 7  j 


r 


t 


3  3 
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xy 
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The  stress  and  moment  resultants  are  related  to  the  reference 
surface  (for  each  region)  strains  and  changes  in  curvature  and  torsion 
by 


where 


(10) 


I  *  1,2, 3, 4 


(1\  -  \-i> 


1=1, 2, 3,4 


(11) 


In  deriving  the  governing  equations.it  is  assumed  that  each  region  of 
the  delaminated  cylindrical  shell  is  symmetric  about  its  mid-surface 


hence  ,  there  is  no  coupling  between  bending  and  stretching  ,  3^- 


0. 


Moreover  ,  only  cases  where  the  fibers  are  in  the  axial  direction  and/or 
in  the  circumferential  direction  are  considered, which  means  that  the 


stiffnesses  A  13  • A  23  ’ Q  13  ,D  23  are  2ero*  T"hus  •  ( 10  ) 

form 
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2.1.3  Equi  I ibr ium  Equations 

The  equilibrium  equations  for  a  thin  cylindrical  shell 
below 


N  +  »  -  0 

x  xy, 

' i,y 


NXV  +  N  ■  0 

xyi,x  yyt 


i.y 


M  +  2  M  +  M  +N/R  +  NW.  + 


xx 


1 , xx  xyi,xy  yyi. 


yy 


yy.  xx.  l.xx 


+  2  N  w,  +■  N  w.  *  0 
xyt  t,xy  yy.  t,yy 


assume  the 

(12-a) 

(12-b) 

are  1 i sted 

(13) 

(14) 


(15) 


2.2  PRIMARY  (MEMBRANE)  STATE  SOLUTION 


As  the  compressive  load  pf  is  applied  (quasi)  statically, a  (membrane) 
primary  state  exists.  Assume  the  initial  curvature  to  be  the  same  for 
all  the  four  regions  i.e.  t  «  R  .  Also, assume  that  a  uniform  radial 
displacement  is  allowed.  Then  ,  the  primary  state  is  characterized 


where 


yp  m  yp  »MP  »NP  >  =>  of 

T/t  xyi  yyi  xyi 


-  H  p 


where  p  is  the  compressive  applied  load  .and  b*»h/t  and  H-H/t 


Furthermore, 


where 


ui  *  -  P,/C  x  ;  vp  -  0 

i  1 


Wi  *  -<Pt/C  )  (A  /A  )  i»i  2  3  a 

*•  1  xx^  xy  yy'i 


X  ■  (a«  -  Wi 


2.3  BUCKLING  EQUATIONS 


The  buckling  equations  are  derived  by  employing  a  perturbation 
technique.  Since  the  perturbed  state  can  be  chosen  to  be 
i nf i ni tess imal 1 y  close  to  the  primary  state  only  ,  first  order  terms  in 
the  admissible  variations  are  retained. 

The  following  expressions  are  substituted  in  the  governing  equations 


a  a  a  pa 

ui  *  ui  +  ai  s  vi  3  vi  ;  wi  3  wi  +  wi 


NP  +  N3  ;  N  *  N3  •  N  *  N3 
xxL  xxL  xy.  xyL  *  ^yy.^  yy. 


Ma  ;  M  -  Ma  ;  M  '  Ma 

**1  i  xyt  yyt  yyt 


Q  -  Q 
Xi  xi 


where  the  super  "a"  parameters  are  referred  to  the  additional 
changes  (that  corespond  to  admissible  variations  uf  ,vf  and  wf  ). 


By  retaining  only  linear 
the  buckling  equations  are: 


terms 


in  the  additional  paramters 


N3  +  N3 

“i.x  xyi,y 


(36 
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xy 


+  N 


t,x  yyt,y 


(37 
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+  2  M 

at,  xy 


+  M 


i ,  xx 


'Tft 


+  /R  -  N 


i.xy  i.yy 


yy 


w. 


xx ^  i,xx 


(38 


In  terms  of  the  in-plane  and  transverse  d i sp I acement  (u, v, w)  ,  th 
buckling  equations  take  the  form 


4  D  wa 
-xxyy  yyt  i.yyyy 


S  “l,—  +  «Va,..,l  Wi,3 

t  Axyt  (\iy-  “*/<»}/»  -  0 


+  8*x  v 


(39 


\  “I  +  ^  ”i,x/R 


(40 


uf  +  L,  vf  =  A  wa  /R 


(41 


where  the  differential  operators  Lpl^and  L  2  ar e  given  by 


L1  *  Axx  +  Asg  b2  /  cy2 

L2  *  (Axy+A3a)  a2  /  axay 

l3  "  A33  a2  /  s*2  +  A  a2  ,  a  2 

yy  '  y 


Equations  (40)  ,(41)  can  also  be  written  as: 


»?  -  a _ C  A...  a3  /a*3  -  Avv.  a3  /a*ay2  3  «!'* 


'i  1  xyt 


yy< 


i.e. 


uj  -  i^tA, ,  U  a3  /a*3  -  a3 


(4 


Li  Vi  "  ^yy/ss^3  /9y3  -(Axyi  +AxyiAssi’AxxlAyyi)a'5  /3x^]  wi^R 


3  «  2, 


i.e. 


•i*  ll‘tcS,*..»3  ^3-(aL_+a, 


yy 


i  i 


xyt  +Axy.Ass."Axx  Ayy.)d3  ^x23y]  w*/Rj 


where 


L. 

x 


(A2  +2a  A  -A  Ayy)^4  /5x2dy2  +A  A  a4  /3y4 
xy  xy  SS  XX  1  i  7  yy  3S.  y 


Lt  I./  -  1 

Equations  (42  )  and  (  43)  are  used  to  eliminate  u^  and  v ^  from  Eq • (  39  ) 
to  obtain  one  (buck  1 i ng)  equation  with  higher  order  derivatives  in  the 
transverse  deflection  w^  . 

Thus, the  buckling  equation  for  each  region  ,  i  ,  becomes 


A  3  /3x  -  (A2  +2A  A  -A  A  )  .S'*  /Sx^y2  +A  A  S4  /Sy4} 

xxt  xy  xy  ss  xx  yy  i.  y  yyi  ssi  y  J 


X  {Dxx  a4  /Sx4  +  2(D  +2Dsa)l  a4  +  d  a4  /ay4  +  n  a2  /ax2>, 


^  /A,,2*,,2  -l.  n  «  a2  2j 


xy  - ss' x  '  yyL 

+  A  A  Ca  -  A2  /A  }.a4  /ax4]  w*/R2  a  0 
ss  xx,.  yy  xy  xx  x,  x 


XX. 
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(44) 


x  x 


i-1, 2,3,4 


The  related  conditions  are  grouped  into  three  categories  :  Boundary 
Conditions  ,  Kinematic  Continuity  Conditions  and  Continuity  in  Moments 


and  Forces. 


Along  each  end  of  the  cylindrical 
have  to  be  satisfied.  For  both 
boundar i es , the  following  possibilities 

a)  Simply  Supported 
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b)  Clamped  Support 
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shell  four  boundary  conditions 
simply  supported  and  clamped 
hold  (Ref  6  )  : 
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where  j*l  at  x*  0  and  j“4  at  x*  l 


2.4  SOLUTION  OF  BUCKLING  EQUATIONS 


The  solution  of  the  buckling  equations  , Eq.  (44  )  .can  be  written  as 


[6]  : 
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Xjj  are  the  characteristic  roots  of  the  buckling  equation,  and  can  be 


found  by  solving  the  following  equation: 
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Note  that  for  continuity  of  displacement  around  the  circumferential 
di rection 
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where  n  is  the  number  of  full  waves  in  the  circumferential  direction. 
Also  note  that  n  must  be  the  same  for  all  regions  ,for  continuity  of 
displacements  to  be  satisfied  . 

Using  Eqs(  42),(  43)  &  (63)  ,  the  displacements  ,  stress  and  moment 
resultants  are  expressed  in  terms  of  ,  S^and  as  follows: 
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The  solution  to  the  buckling  equations,  Eqs  (63) .  requires  knowledge 
of  32  constant  (Aij, i*1 ,2,3*4  , j«l ,2, 3. • • *8) .  There  exist  32  boundary 
and  continuity  conditions,  which  are  homogeneous  in  u^  ,  and  w^  and 
their  space  derivatives.  These  consist  of  eight  boundary  conditions, 
four  at  each  end,  Eqs  (45)  or  Eqs  (46) ,  sixteen  kinematic  continuity 
conditions,  Eqs  (47) - (54) ,  and  eight  continuity  conditions  in  moments  and 
forces,  Eqs  (55)  ~  (62)  . 

The  use  of  the  above  system  of  boundary  and  continuity  conditions 
yields  a  system  of  linear,  homogeneous,  algebraic  equations  in  A^j  . 
The  resulting  system  of  equations  takes  the  form 

[C]  C  x3-  [0]  (71) 

where  [C]  is  a  32x32  matrix  and  { x}  an  1x32  (column)  matrix.  The 
elements  of  matrix  C  are  given  in  appendix  A.  For  a  nontrivial  solution 
to  exist,  the  determinant  of  the  coefficients  must  vanish.  The 
determinant  contains  geometric  parameters  (L,R,a,h),  material  parameters 
(Eh  ,£22  »<»i2  »vi2^  »  the  aPP''ec*  load,  p,  and  the  circumferential  wave 
number,  n. 

For  a  given  geometry  and  material  properties,  the  circumferential 
wave  number  (n)  is  varied  and  the  corresponding  load  which  makes  the 
determinant  vanish  is  obtained.  The  lowest  of  these  loads  is  the 
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2.5  RESULTS  AND  DISCUSSION 

Buckling  analysis  of  delaminated  cylindrical  shells  involves  many 
geometric  and  material  parameters.  Also,  as  is  well  known,  the 
different  inplane  boundary  conditions  have  a  significant  effect  on  the 
critical  load  of  the  cylindrical  shells.  Therefore,  general  results, 
which  account  for  all  the  effect  of  the  different  parameters  on  the 
critical  load  of  the  delaminated  cylindrical  shell,  require  several 
extended  studies.  Accordingly,  the  results  presented  herein  are 
restricted  to  specific  numerical  examples  in  order  to  illustrate  the 
applicability  of  the  present  model. 

In  generating  results  for  the  present  model,  it  is  assumed  that 
there  is  no  contact  between  the  two  regions  across  the  delamination  at 
the  instant  of  buckling.  Otherwise  the  problem  has  to  be  treated  as  if 
a  cylindrical  shell  is  resting  on  an  elastic  foundation  (not  attached  to 
the  foundation).  A  model  for  accounting  for  thin  behavior  is  suggested 
in  Appendix  B. 

For  a  given  configuration,  the  critical  load,  p,  is  determined  as 
the  smallest  eigenvalue  that  makes  the  determinant  of  the  coefficients 
of  the  system  of  homogeneous  algebraic  equations  (  see  Eq . 7 1 )  equal  to 
zero.  The  circumferential  wave  number  is  varied  in  order  to  determine 
the  smallest  of  the  eignvalues. 

Results  are  generated  for  a  cylindrical  shell  made  up  of  isotropic 
laminate,  and  for  both  simply  supported  and  clamped  boundary  conditions. 
Since  for  each  type  of  supports  (simple  and  clamped  supports)  there  are 
four  types  of  boundary  conditions,  results  are  generated  for  the  weakest 
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and  for  the  strongest  configurations,  which  are  known  as  SSI  (w  *  w^^  * 
axx  *  TXy  ■  0)  and  CC4  (w  *  w  >x«  u  ■  v  ■  0)  ,  respectively.  The 
boundary  conditions  are  assumed  to  be  the  same  at  both  ends,  i.e.  both 
ends  are  SSI  or  CC4.  The  dimenstions  of  the  cylindrical  shell  are  such 
that  L/R*5i  and  R/t*30>  where  L,  R  and  t  are  the  length,  radius  and 
thickness  of  the  shel 1,  respectively. 

Figure  3  shows  the  effect  of  delaminatiand  length  and  2-position  on 
the  critical  loads  of  a  simply  supported  cylindrical  shell  (SSI).  The 
delamination  is  assumed  to  be  located  symmetrically  with  respect  to  both 
ends  of  the  shell.  The  critical  loads  are  normal i2ed  with  respect  to 
the  critical  load  of  the  classical  theory.  The  figure  shows  that  for  a 
delamination  in  the  middle  surface  of  the  cylindrical  shell,  the 
delamination  has  a  negligible  effect  on  the  critical  loads  as  long  as 
the  delamination  is  far  from  the  edges  (  effect  of  position  with  respect 
to  the  edges  will  be  discussed  in  another  figure  ).  As  the  delamination 
moves  away  from  the  mi dd 1 e-  surface  of  the  shell,  its  presence  becomes 
important.  For  a  delamiation  thickness  h  ■  0.3,  the  delamination  has  no 
effect  on  the  critical  loads,  as  long  as  the  delamination  length  is 
small  (a  <.04) .  As  the  delamination  length  increases  the  critical  load 
decreases  and  it  approaches  asymptotically  a  value,  which  is, 
approximatly  about  60%  of  the  critical  load  of  the  perfect 
configuration.  As  the  delamination  becomes  thinner  and  thinner,  a  sharp 
drop  in  the  critical  load  is  noticed  at  a  very  small  delamination  length 
(  a  -  0.01  ) .  This  drop  in  the  critical  load  continues  till  it  reaches 
a  constant  value  (  20%  of  that  of  the  perfect  shell  )  ,  at  a  value  for 
the  delamination  length  paramter  of  I  ■  0.1  . 


The  effect  of  lengthwise  delamination  position  on  the  critical  loads 
is  also  studied.  Figure  4  presents  results  for  1  -  0.1  .  For 
delamination  thickness  h  $  0.3.  the  delamination  position,  l,  has  no 
appreciable  effect  at  on  the  critical  loads.  For  thicker  delaminations, 
h  >  0.3,  the  position  of  delamination  has  an  important  effect  on  the 
critical  load,  and  this  effect  increases  as  the  delamination  thickness 
increases.  For  instance,  a  reduction  of  about  one  third  in  the  critical 
load  is  noticed  for  a  delamination  near  the  edge  (by  comparison  to  the 
critical  load  of  a  delamination  far  from  the  edge  ;  ^  >  .08  ),  for  a 
delamination  thickness  h  *  0.5  •  The  position  effect  decreases  as  the 
delamination  thickness  decreases.  Figure  5  shows  the  effect  of  both 
delamination  length  and  position  on  the  critical  load  of  a  delaminated 
cylindrical  shell,  with  the  delamination  located  at  the  mid  surface  of 
the  shell  (h  ■  0.5).  For  a  moderate  delamination  length,  0.1<a  <0.2,  it 
seems  that  the  effect  of  delamination  position  is  approximately  the  same 
for  bothe  values  of  delamination  lengths  (  a  *  0.1  and  0.2  ).  For 
shorter  delamination  length,  the  critical  load  decreases  as  the 
delamination  moves  away  from  the  ends  of  the  shell  and  it  reaches  a 
minimum  at  a  value  of  l  in  the  range  0.04  >  l  >  0.02.  The  exact  value 
of  ^  at  which  the  minimum  is  reached  depends  on  the  delamination  length 
parameter,  a.  As  the  delamination  length  parameter  becomes  shorter  and 
shorter,  (a  ^  0.02),  the  effect  of  delamination  position  can  be  ignored. 
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A  similar  set  of  curves  are  obtained  for  clamped  boundary 
conditions.  The  effect  of  delamination  length  on  the  critical  loads  of 
a  symmetrically  delaminated  cylindrical  shell,  for  different  values  of 
the  delamination  thickness,  is  shown  in  Figure  6.  This  figure, 
illustrates  that,  only  for  very  small  values  for  the  del  ami n?r i un  length 
(a  «  .01),  the  delamination  has  no  significant  effect  on  the  critical 
loads,  regardless  of  the  value  for  h.  For  larger  values  for  the 
delamination  length  parameter,  a,  a  noticeable  drop  in  the  critical 
loads  is  observed.  For  example,  a  drop  of  about  one  half  of  the 
critical  load  (  of  the  perfect  shell  )  is  observed  at  a  delamination 
thickness  parameter,  h  ■  0.5.  The  drop  of  the  value  of  the  critical 
load  is  more  severe  for  thinner  del  ami nat i ons .  For  instance,  the 
critical  load  is  only  10%  of  the  critical  load  of  the  perfect  shell  for 
delamination  thickness  h  *  0.1  with  delamination  length  a>.07. 

Unlike  the  simply  supported  case,  the  position  of  delamination 
w.r.t.  the  edges  of  the  shell  with  clamped  ends,  has  virtually  no 
effect  on  the  critical  loads,  for  delamination  length  a  *  0.1, 
irrespective  of  its  thickness,  as  shown  on  Figure  7* 

The  effect  of  delamination  position  on  the  critical  loads  of  a 
delaminated  cylindrical  shell  with  the  delamination  located  at  the 
middle  surface  of  the  shell,  for  various  values  of  delamination  length 
is  presented  on  Figure  8.  It  is  observed  that  for  moderate  delamination 
lengths  (a  *  0.1  and  0.2)  the  position  of  the  delamination  has  no  effect 
on  the  critical  loads  of  the  shell.  On  the  other  hand,  for  short 
delaminations  (a  ■  .05),  the  critical  load  decreases  as  the  delamination 
moves  away  from  the  edge  of  the  shell. 
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A  comparison  of  the  results  on  Figures  3  and  6,  shows  that,  for 
moderate  and  large  delamination  lengths  (a  >  0.1),  a  delaminated 

cylindrical  shell  has  the  same  critical  load,  irrespective  of  the 
boundary  conditions  (simply  supported  or  clamped)  and  this  critical  load 
paramter  is  approximatly  equa 1  to  H.  The  explanation  of  such  results  is 
that  for  such  delamination  length  paramter  (a  >  0.1),  the  thin  region 
(delaminated  part)  buckles  while  the  main  body  of  the  cylindrical  shell 
remains  straight  (has  no  appreciable  bending  deformation),  i .e.  the 
thin  part  buckles  as  if  it  had  a  clamped  boundary  condition, 
irrespective  of  the  type  of  supports  of  the  main  cylindrical  shell. 

Figure  9  shows  a  plot  of  the  critical  load  parameter  for  region  3 
(the  thin  region)  versus  the  characteristic  length  of  region  3  (L3  *  a  / 

Jr  h  )  for  a  very  thin  and  long  delamination  (h  •  .02,  a  ■  0.6),  for 
clamped  boundary  conditions.  The  main  purpose  of  presenting  this  curve 
for  this  configuration  (very  thin  and  very  long)  is  to  check  a  limitting 
case,  which  is  similar  to  the  thin  film  analysis  of  the  delaminated 
plate.  As  expected  the  obtained  curve  is  similar  to  the  one  obtained  by 
Hoff  [6]  for  the  critical  loads  of  a  cylindrical  shell  with  CCU  boundary 
cond i t i ons . 
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2.6  CONCLUDING  REMARKS 


The  investigation  of  the  problem  of  buckling  of  delaminated 
cylindrical  shells  leads  to  the  following  observations: 

1)  For  symmetric  delamination  (delamination  located  symmetrically  w.r.t. 
both  ends  of  the  shell)  the  critical  load  is  independent  of  the  boundary 
conditions  as  long  as  the  delamination  length  parameter,  a,  is  greater 
than  0.1  .  For  this  case,  1  >  0.1,  the  critical  load  depends  only  on 
the  delamination  thickness  parameter,  h. 

2)  For  clamped  boundaries  (CC-U)  ,  the  position  of  delamination  w.r.t. 
the  ends  of  the  cylindrical  shell  has  no  effect  on  the  critical  load  for 
moderate  and  large  values  of  the  delamination  length  parameter,  I>0.1  . 
As  the  delamination  length  becomes  smaller,  the  effect  of  position  of 
delamination  becomes  more  pronounced,  and  the  critical  load  increases  as 
the  del  ami  nation  moves  towards  the  clamped  ends  . 

3)  For  simply  supported  boundaries  (SS-1)  ,  the  position  of  delamination 
w.r.t.  the  ends  of  the  cylindrical  shell  has  a  significant  effect  on 
the  critical  load,  even  for  small  values  of  the  delamination  length 
parameter.  Unlike  the  case  of  clamped  boundaries,  for  simply  supported 
boundaries  (  with  a  *  0.1  ),  the  critical  load  decreases  as  the 
delamination  moves  towards  the  edges  of  the  shell  boundary  as  long  as 
the  delamination  thickness,  h,  is  greater  than  0.3.  otherwise,  the 
change  of  delamination  position  has  a  negligible  effect  on  the  critical 
load  for  smaller  values  of  the  delamination  thickness,  h  <  0.3  • 


APPENDIX  A 


CHARACTER  I  ST  1C  EQUATIONS  OF  DELAMINATED 
CYLINDRICAL  SHELLS 


The  solution  to  the  buckling  equation  ,  Eq. (44  ),  can  be  written  in 


the  form  : 
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J-1  J 
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Use  of  the  proper  boundary  and  continuity  conditions  leads  to  the  system 
of  linear  ,  homogeneous  algebraic  equations 

[c]{x]  -  0 

where  [cl  is  a  32x32  matrix  and  £ x]  a  1x32  (column)  matrix, 
r  T 

[x]  -  [  \i’Ai2’Al3,Alb’Al5,Al6,Al7,Aia,A2l,A22,A23’A2U,A2S’A26'A27’A28’ 

A31 ’ A32 ’ A33 ’ A34 ’ A35 ’ A36 ’ A37 ’ A38 ’ A41 ’ A42 ’ A43 ’ A44 ' A45 ’ A46 ’ A47 ’ A48^ 


The  general  element  of  the  matrix  [C]  is  c 

ij 


In  defining  the 


elements  two  parameters  (r,q)  are  used,  and  defined  as 

j  a  1,2,3  ...  8  ,  r*l  and  q=j 

j  -  9,10,11,..  ,  r»2  and  q-j-9 

j  *  17, 18,..,  24  ,  r-3  and  q=»  j -16 

j  *  25, 26,... 32  ,  r-4  and  q-j-24 


Also  we  have. 


ctj-  0  i-1,2,3,...  16  and  j  -25,26, 

and  for  i-17,18,19. . .32  and  j  -1,2,3, 
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Moreover,  the  following  parameters  are  defined  as 
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The 

elements  of  Che 

matrix  [ c]  are  obtained  by  multiplying  each 

of  Che 
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The  above  elements  (c^)  corespond  to  the  clamped  supports  case  of 
CC-4  .  For  simply  supported  boundaries  of  SS-1  the  elements  with  the 
super  star  must  change  to  , 
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APPENDIX  B 


COLUMNS  PARTIALLY  SUPPORTED 


BY  ELASTIC  FOUNDATION 


Structural  elements  supported  elastically  along  part  or  along  their 
entire  lengths  are  used  very  offen  in  structure  conf i nurat i ons .  The 
elastic  supports  may  be  a  foundation  as  in  railroad  tracks,  or  adjacent 
elastic  structural  elements  such  as  in  stiffened  plates  and  shells  [7]. 
Moreover,  elastic  supports  can  also  be  found  in  delaminated  structures 
when  the  delaminated  region  (or  part  of  it)  deforms  towards  the  main 
body  of  the  shel  1  , 

Although  there  are  many  models  to  represent  the  elastic  support 
(which  will  be  referred  to  as  the  elastic  foundation  from  here  on), 
[8-13],  the  simplest  of  all  these  models  is  known  as  the  Winkler 
foundation  [8]. 

The  considered  elastic  foundation  is  assumed  to  be  of  the  Winkler 
type  [8],  where  the  foundation  is  approximated  by  a  series  of  closely 
packed  linear  springs,  and  the  foundation  reaction  at  any  point  is 
proportional  to  the  deflection  at  that  point. 

A  column  of  total  length  L,  has  parts  of  length  /  ,  which  are 
surrounded  by  the  elastic  foundatiopn,  so  that  the  elastic  foundation 
has  the  same  lateral  deflection  as  the  column,  for  these  parts.  The 
rest  of  the  column  is  not  supported  elastically  (see  Figure  B-l)  . 


The  buckling  equations  that  govern  the  behavior  of  the  two  regions  (with 
and  without  the  elastic  foundation)  of  the  axial  1 y  column  are  : 


p  a 

w  -t w  =  0 

,xxxx  El  ,xx  El 


v  +  -£•_  w  =0 

,xxxx  El  ,  xx 


0  <  x  <  2 


(B-l ) 


(B-2) 


2  <  x  <  L 


wher  9  is  Che  modulus  o£  elastic  foundation 


The  solution  of  the  buckling  equation,  Eq.(B-l),  takes  one  of  the 
following  three  forms  .depending  on  the  value  of  the  load  parameter  (7 
■  p/2 \/^TT)  ,  (for  details  see  Ref. 39) 


Case  I  : 


Y  >  1 

*  bj  cos  k^x  +  b2  sin  k^x  +  by  co3  k2x  +  b^  sin  k2x  (3-3i) 


where 


Case  II  : 


k,  =•  (0/EI)*  cy  -  Jy2- 1)* 


y  »  l 


Wi  =  b,  cos  k.x  +  b„  sin  k„x  +  b,  x  cos  k.x  +  b,  x  sin  k.x 

*•1  3  2  3  3  34  3 


(B-3 ii) 


where 


k3  =»  (fl/EI)' 


Case  III  : 


Y  <  1 


,  (s+Ltlx  ,  , 

u.  zt  h  s  x  h  a 

1  ui  2 


where 


(s+it)x+  b  e(-s+it)x  +  b  e(s-it)x 
3  4 


s-O/EI)  v'(l-Y)/2 
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The  solution  of  the  buckling  equation,  Eq.  (B-2)  ,  is  of  the  form 


w* 


cos  kx  +  3^n  kx  +  x  +  34 


L  >  x  >  l 


where 


k  =*  J  p/EI 


The  related  boundary  and  continuity  conditions  are 
Boundary  conditions 

Either  Or 


(B-4) 


w  =*  0 


w  *  0 

.x 


Continuity  Conditions  ^at 


w  + 

,  XXX 


p_ 

El 


w 

,x 


0 


(B-5) 

w  *  0 
,xx 


Wl,xa  W2,X 

wl,xx=  W2,xx 

w,  =  w_ 

1 , xxx  2 , xxx 


(B-6) 


The  solution  of  the  buckling  equations,  Eqs.(8-3)  and  (B-4-)  requires 
knowledge  of  eight  constants  (b(*  ,a-  ,  i*l,2,3.4).  There  exist  eight 
boundary  and  continuity  conditions  .  The  use  of  the  boundary  and 
continuity  conditions  yields  a  system  of  linear  ,  homogeneous  , 
algebraic  equations  in  bj  and  a-  .  For  a  nontrivial  solution  to  exist 
the  determinant  of  the  coefficients  must  vanish. 
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Results  are  obtained  for  axially  compresed  columns  with  simply 
supported  boundary  conditions.  The  same  solution  scheme  is  also  used  to 
obtain  critical  loads  for  a  simply  supported  column  with  one  (two)  third 
of  its  total  length  elastically  supported.  Moreover,  results  are 
presented  for  columns  divided  into  four  equal  parts,  two  of  them  being 
elastically  supported. 

The  obtained  results  ar^  presented  graphically  on  Figures  (B-2)  and 

(B-3).  On  Figure  (B-2)  the  critical  load,  normalized  w.r.t.  p  (p  ■ 

E  E 

tt2ei/l2  ),  for  a  column  with  part  of  length  l  ,  attached  to  an  elastic 
foundation,  while  the  rest  of  the  column  is  not.  As  the  length  of  the 
portion  of  the  column  which  is  attached  to  the  elastic  foundation 
increases,  the  critical  load  increases.  Of  course  all  the  curves  lie 
between  the  two  limiting  cases,  the  straight  line,  P  ■  1.0,  which 
corresponds  to  the  case  where  there  is  no  elastic  foundation  at  all,  and 
the  curve  presented  by  the  broken  line  .which  corresponds  to  the  case 
where  the  column  is  fully  supported  by  the  elastic  foundation. 

Figure (B-3)  presents  the  critical  load  parameter,  Pcr  ,  versus  the 
normalized  modulus  of  foundation,  f  .  On  Figure (B-3),  the  column  is 
assumed  to  be  devided  into  equal  parts  (two, three  and  four)  and  one 
(two)  part(s)  of  the  column  is  elastically  supported.  These  curves  are 
obtained  primarily  to  have  an  (approximate)  idea  about  what  changes  are 
expected  to  happen  to  the  critical  load  of  the  delaminated  cylindrical 
shell,  if  there  is  a  contact  between  the  two  regions  across  the 
delamination.  As  the  results  of  Figure(B-3)  show,  the  value  of  the 
critical  load  depends  on  the  wave  number,  the  dimension  of  the 
structural  element,  as  well  as  on  the  (equivalent)  modulus  of  foundation 
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1 .  INTRODUCTION 


The  stability  of  multi-annular  plates  is  being  considered.  Most 
researchers  Investigated  the  buckling  of  a  single  annular  plate  (with  or 
without  a  rigid  inclusion)  when  subjected  to,  primarily,  a  uniform  stress 
field. 

The  elastic  stability  of  an  annular  plate  subjected  to  shear  forces, 
distributed  along  the  edges,  was  considered  by  Dean  [1].  Stability  due  to 
bending  moment  caused  by  initial  stress  was  considered  by  Willers  [2]. 
Federhofer  C 3 , ^ D ,  Olsson  [5],  and  Eggar  [6]  dealt  with  the  buckling  of 
annular  plates  of  varying  thickness.  Stability  analyses  under  a  uniform 
compressive  load  were  considered  by  Olsson  [7],  Schubert  [8]  for  clamped 
and  simple-supported  plates,  respectively,  and  by  Yamaki  [9]  for  various 
boundary  conditions.  Yamaki  showed  that  'under  some  B.C.'s  the  axisymmetric 
mode  does  not  yield  the  lowest  critical  load,  and  therefore  higher  modes 
must  be  investigated. 

Meissner  [10]  dealt  with  axisymmetric  buckling  of  an  annular  plate, 
3implys>3upported  or  clamped  at  the  outer  boundary  and  free  at  the  inner  one. 
Most  of  the  above  references  used  Bessel  functions  to  solve  the  buckling 
equation.  Mansfield  [11]  dealt  with  an  infinite  circular  plate  with  a 
hole  supported  along  two  concentric  circles  and  subjected  to  radial  inplane 
loading  along  the  inner  circles.  The  solution  achieved  was  of  the  closed 
form.  The  same  type  of  problem  was  considered  by  Bareeva  [12-14]  and  by 
Lizarev  [15-17].  Bareeva  [12-14]  analyzed  a  single  annular  plate  with  a 
very  rigid  ring  at  the  outer  boundary.  A  power  series  solution  was  used. 
The  results  are  given  for  the  first  and  second  modes  only.  Lizarev  [  1 6 **  1  7 ] 
employed  the  hypergeometric  function  for  formulating  the  general  solution. 
Results  are  shown  only  for  the  first  and  second  modes.  A  comparison  of 


theoretical  and  test  results  was  presented  by  Majumdar  [18]  for  a  single 
annular  plate.  The  solution  for  the  first  two  modes  was  achieved  by  using 
Bessel  functions.  For  the  other  modes,  a  Galerkin  procedure  was  employed. 
The  convergence  of  the  Galerkin  solution  Is  limited  to  0.2  <  r  inner f 
r  outside  ^  0.6. 

Tie  results  appearing  In  the  above  references  are  usually  limited  to 
the  first  two  modes  and  are  only  applicable  to  a  single  annular  plate.  The 
results  are  so  limited  because  of  difficulties  in  the  mathematical 
formulation  and  the  employed  solution  schemes. 

The  present  investigation  deals  with  the  buckling  problem  of  a 
mult  1 ►annular  circular  plate  with  different  material  properties  and 
thicknesses.  The  in**plane  loading  is  axisymmetrlc  and  many  combinations  of 
transverse  and  inplane  boundary  conditions  are  considered. 

The  stability  analysis  employs  a  power  series  type  of  solution  to 
describe  the  buckling  modes  (functions).  The  power  series  method  Is  used 
because  of  its  applicability  to  the  various  B.C.'s  and  different  properties 
and  because  it  enables  a  general  formulation  for  the  multi-annular  problem 
(buckling  analysis). 

The  described  formulation  and  solution  procedure  are  general.  Results 
are  presented  only  for  solid  and  annular  circular  plates.  The  manuscript 
outlines  the  basic  assumptions,  the  solution  scheme,  and  the  principal 
parameters  affecting  both  the  critical  load  and  the  type  of  solution.  The 
complete  mathematical  formulation  is  decrlbed,  with  detail. 

The  mathematical  formulation  Is  based  on  the  following  simplifying 
assumptions: 

(I)  The  material  is  linearly  elastic,  satisfying  Hooke's  Law. 

(II)  The  properties  are  constant  within  each  annular  section. 


(Hi)  The  stresses  are  smaller  than  the  proportional  limit. 

(iv)  The  usual  Klrchhoff«Love  hypotheses  are  applicable. 

(v)  Each  annular  section  of  the  plate  is  thin  (the  ratio  of  thickness  to 
outer  radius  is  small  by  comparison  to  one). 

(vl)  The  deformation  gradients  and  the  rotations  are  small. 

2.  ANALYSIS  OF  A  MULTI -ANNULAR  PLATE 
The  buckling  analysis  of  multi  annular,  thin  plate,  subjected  to 
axisymmetric  loading,  is  described  by  the  solution  of  the  following 
equations  for  each  annular  part  (see  Appendix  A. 3). 


(rN°  )  -  N°  -  0 

rr  ,r  96 


0) 


D  v\f 


N99,9  =° 

N°  w,  ♦  N®  (  ^  ) 

rr  ’rr  99  v  2  r  ' 


(2) 


where  N^,  N°Q  are  primary  state  stress  resultants  in  the  radial  and 
circumferential  directions,  respectively;  w  is  the  transverse 
displacement;  (  ),  ,  (  ),  are  partial  derivatives  with  respect  to  radial 

or  angular  coordinates;  and  D  is  the  flexural  rigidity. 

Eqs.  (1)  govern  the  primary  state.  The  other  equation,  Eq.  (2),  is 
the  buckling  equation  and  is  derived  through  linearization  of  the 
transverse  equilibrium  equation.  The  primary  state  is  first  solved  for  the 


entire  plate.  Then,  use  of  these  solutions  (N°r,  N°0),  in  Eq.  (2),  leads 


to  the  buckling  analysis. 


2.1  Primary  State 

The  governing  equilibrium  equations  of  the  primary  state,  for  each 
annular  section,  are  given  by  Eqs.  (1).  The  in^plane  stress  distribution 
(see  Fig.  1)  is  determined  through  inf-plane  equilibrium,  Eq.  (1),  and  the 
following  in*plane  displacement  and  force  continuity  conditions  of  the 
adjoining  boundaries  (of  neighbouring  annular  plates): 


Jii  “  ai+1  i 


N°  =  N° 
rril  rri+1 i 


where: 


,Jit*  '-Ui  +  1  are  the  radial  displacements  of  plates  i  and  i  +  1  at  the 
common  joint,  i. 

N  °  .  N  °  are  the  radial  stress  resultants  at  joint  i  of  plate  l  and  i 
ii  rri+1i 


The  stress  distribution  is  derived  by  using  the  stiffness  approach 
(for  details  see  Ref.  20).  The  equations,  using  matrix  notation,  are: 


S  u  -  p 
- u-e 

where:  S  is  the  stiffness  matrix  [s.,]  £,j  =  0,  1...,  N-1 ; 

tj 

u  is  the  displacement  vector  at  the  various  joint,  [u^]; 

£e  is  the  boundary  stress  resultant  vector  at  various  bounding 


circles  [p  ];  and  i,j  are  integer  subscripts 
el 


The  sign  convention  for  displacements  and  forces  is  described  on  Fig.  1. 
The  stiffness  matrix  is  not  symmetric  and  the  relation  between  the 
off-diagonal  terms  is: 
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*  x.j  J  jx 


where  r^  are  the  radii  at  joints  l  and  j 


The  stress  field,  for  an  annular 


plate  i,  is  given  by  (see  Ref.  20) 


N°  =  -N  +  N  /r2 
rri  °i  °Ci 


80.  »  -N  -  N  /r 
i  ot  oGi 


N  Q  -  0 
r0i 


N  *  (N  B?  ^  N  )/  ( ef-1 ) 

°i  rrii-1  1  rrii  1 


Noc  *  S?ri  (Nrr  "  Nrr  )/(Bi  "  15 
i  ii-1  ii 


(6) 


where:  N°  ,  N  °  ,  N  °  are  the  radial,  circumferential  and  shear 
rrt  90  A  rOj 

stresses  in  plate  i. 

6^  -  ri-i/ri  “  ratio  of  outer  to  inner  radii  for  plate  i. 
r'i’i1  ’  ri  *"  outer  and  *nner  radii  of  plate  i. 

The  stress  field  in  the  central  part  is  uniform  and  equals: 


rr. 


08. 


-  N 


rr 


N-1  N-1 


(7) 


n  4 


where  N  is  positive  if  acting  outward. 

rVl  N-1 


•  )  Forc«  b>  Deformations  and  Beomtry 

Fig.  1  Forces  and  Deformation  in  Annular  Plate 

2.2  Buckling  State 

The  buckling  equation  for  each  annular  or  circular  part  is  given  by 
Eq.  (2).  The  detailed  solution  of  the  equation  using  a  power  series 
procedure  is  described  in  Chapter  3.  The  characteristic  equation  is 
obtained  through  the  satisfaction  of  the  boundary  conditions,  continuity  at 
common  joints  and  kinematic  constraints,  if  they  exist  (see  Fig.  2). 
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Fig.  2  Loading  and  Section  of  the  plate 


The  continuity  conditions  which  must  be  satisfied  at  joint  i,  r  =  are: 


w  *  w 
ii  i+li 


W  3  w 

ii,r  i+1 i  ,r 


M  »  M 
ii  ii  +  1 


Qii  =  _Qi  i+1 


(8) 


where  (see  Ref.  19) 
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rr 
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and  D4  -  Eir'i/12  (1  “  v^)  ; 


(  )^  (  is  a  "plate  i+1"  parameter  evaluated  at  joint  i, 


(The  first  index  designates  the  plate  and  the  second  the  joint); 


Mj»r,  Qr  are  the  bending  moment  and  shear  force  resultants 
jpectively,  and  N°  is  the  radial  stress  resultant  of  plate  i  at  joint 


The  conditions  of  the  boundaries  are  as  follows: 
Edge  boundaries: 

Clamped:  w  =  w>r  *  0 

Simply-’supported:  w  =  M  =  0  [see  Eq.  (8))] 


Free: 


M  =  Q  »  0  [See  Eqs.  (8)  -  (10)] 


Free  to  deflect  not  to  rotate:  w,^  «  Q  =*  0  [see  Eqs.  (8)  and  (10)]. 

Free  to  rotate  but  elastically  supported  in  transverse  direction: 

_  q-i<  w  =  0  (k  is  the  extensional  spring  constant  in  the  transverse 

w  w 

direction). 

Elastically  restrained  against  rotation  and  transverse  translation: 

"  “  krot“V  ‘  Q  "  R-  “  '  ° 

(krot  is  the  rotational  spring  constant). 

If  any  kinematic  constraints  exist  at  an  Interior  joint,  i,  the 
following  equations  replace  Eqs.  (8). 

Transverse  (knife  edge)  support:  w^  *  wi+11  =  Wii,r  =  wi+11fr) 

Mtt  =  Mii+1 
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Completely  fixed  : 
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,  _  M  _  M 

1+1  1  * 

Is  the  extensional  spring  constant  at  joint 


Transverse  extensional  elastic  support:  ^  »  wi+1i  r* 
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1  In  the  transverse  direction). 


The  deflection  of  each  annular  or  circular  part,  w^ ,  Is  described  by: 


wl(r,9)  -  I  [  B)nlF,nl<r>  ♦  B2nlF2nl(r)  *  B  F  Cr) 
n«0 


+  B4nlF4ni(r)]  003  n0  (11) 

where: 

Fkni(r)  denote  the  four  independent  solutions  for  the  plate,  corresponding 
to  the  wave  numbers,  n;  k  assumes  the  values  of  1,2,3  and  4  to  the  four 
independent  solutions.  Bknl  denote  constants  to  be  determined  through 
boundary  and  joint  conditions. 

The  characteristic  equation  is  solved  by  using  the  Newton-Raphson 
procedure.  In  cases  where  the  convergence  is  poor  a  Regula^Fals i  or 
bi-section  method  is  used  (See  [21]). 

The  details  and  a  clarification  for  the  employed  symbols  Is  presented 
In  the  subsequent  sections. 

3.  SOLUTION  OF  THE  BUCKLING  EQUATIONS 
A  power  series  solution  is  used  to  solve  the  fourth  order  ordinary 
differential  equation,  derived  after  separation  of  variables.  This  type  of 
solution  Is  chosen,  because  of  the  versatility  of  the  solution,  and  its 


applicability  to  the  various  boundary  and  joint  conditions,  and  the  various 


types  of  loading.  The  separation  of  variables  is  achieved  by  assuming  the 
deflection  to  be  of  the  form, 


w^r,©)  »  l  FnJ(r)  003  n0 
n-0 

where  r,  8  are  the  radial  and  angular  coordinates,  and  F  ^(r)  are  the 
deflection  functions  in  the  radial  direction  (corresponding  to  plate  i  and 

wave  number  n). 


4 

Substitution  of  Eq.  (6)  into  Eq.  (2),  multiplication  by  r  and  omittng  the 
trigonometric  functions,  yields, 
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where  (  )'  -  d/dr;  af*  No./D.  (see  Eq.  [6]);  and  8  =■  N  /N 

i  i  i  c^  oc^  o^ 

A  power  series  solution,  of  the  form  given  below,  is  assumed 


F  (r)  -  l  A  ..  rJ  +  3ni 
ni  nij 


(14) 


Substitution  into  Eq.  (13)  yields 
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This  equation  when  expanded  yields  terms  multiplying  powers  of  r  of 


3ni+j 

the  form  r  ,  with  j  -  0,1, ...2.  The  multiplying  terms  contain  Ani0 

for  j  -  0;  An|i  for  j  =  1;  and  Anjj  and  Anjj*.2  for  all  other  j-values 

(2,3,.*.).  It  can  easily  be  shown  that  one  solution  is  obtained  by 

requiring  An^o  *  0  and  Ani1  *  0«  and  another  solution  is  obtained  by 

requiring  Anjg  *  0  and  An£-j  *  0.  Both  lead  to  the  same  solution  for  the 

buckling  problem. 

Thus,  for  j  »  0  the  indicial  equation  is: 

(s  .-1  )**  (2  +  2n2  +  B  .a  2)(s  «1)2+  [(n*-1)2+  B  a2(1~n2)]  -  0  Cl 6) 

ni  ci  i  ni  c^  i 

Note  that  this  equation,  Eq.  (16),  and  all  those  corresponding  to  other 

values  of  j  (*  0)  are  quartic  algebraic  equations  in  sn<[. 

The  four  roots  of  this  particular  quartic  equation,  Eq.  (16),  are: 


o*"2  * 


B=i“i 


)  ♦ 


I""2*  *  -TT- 


Y  +  5  . 
ni  ni 


<snl-,>l.4-  <'*n2  * 


eoA 


.2  4 

„_2„  2  6ciai 

2n  *  __ 


Y  5  . 
ni  ni 


(17) 


where  Y^  13  equal  to  the  terms  in  parentheses,  and  6^  is  equal  to  the 
terms  under  the  radical  sign  of  the  right-hand  side  of  Eqs.  (17).  The 
recurrence  formula,  in  the  general  solution  [for  all  k  -  see  (Eq.  11)]  is: 


6? 


alC(j+3kni^2  "  n2]  A 


knlje2 


'<nlJ  [U*W>V  s<2  *  2"2*  eW)  W*knl*  1)2(aknr1)2:i 


(j  -  2,4...  even) 

(k  -  1,2, 3, 4) 


(18) 


The  various  types  of  solution  depend  on  the  type  of  roots  (real,  equal  or 
complex).  The  general  solution  (for  the  case  of  four  distinct  roots); 


Fkni(r) 


l 

j-0 


knlj 


J  +  3 


knl 


(k  »  1,2, 3, 4) 


(19) 


where  n  is  the  wave  number; 

i  refers  to  annular  plate  1; 
k  identifies  the  four  independent  solutions; 
j  refers  to  the  j'-term;  and 

3kni  denote3  the  kr'h  root  of  the  indicial  equations,  Eq.  (16),  see 
Eqs.  (17),  for  plate  i  and  wave  number  n. 


3.1.  Types  of  Roots,  s^n; 

The  roots  of  the  indicial  equation,  Eq.  (16),  can  be  either  real, 
equal  or  complex  conjugate.  First,  we'll  determine  the  various  types  of 
the  roots,  and  later  we’ll  discuss  the  intervals  of  validity.  In  general, 


ee 


the  roots  can  be  written  as 


(20a) 


where 


3.  .  ■  sf  .  +  1  3^  . 

knl  knl  knl 


Re  Csnl] 


1  *  P 


lm  C3nl]  '  4l1 


'  *  P 


(20b) 


1  +  P 


*  4 

1  "  P, 


R  I  R  X 

The  values  of  Pnl1 f  Pnl1  ,  Pnl2»  pni2  dePend  on  the  values  of  6ni  [see  Eqs. 
(17)].  For  6  >  0  the  quantities  In  Eq.  (20)  are  as  follows:  First,  let 


+  « *4  •  <0. 


Knl1  -  'nl  •  uni  ‘  *ni2  "  Ynl  “  6nl 


If  *nll  *  0,  then  p"n 


j^nil  : 


pnir  0 


if  ^ni1  <  0,  then  pni1  -  0 


if  *nl2  2  0,  then  pni2 


:  pni2  “  0 


If  *nl2  <  0,  then  pni2  -  0 


Pni2  “  pnlT 


Moreover,  for  5nj  -  0  we  have 


if  Y  .  >0,  then  p 
ni  Kni1 


f \Tl  ■■  Pni, 


If  Ynl  -  0.  then  pnl, 


pni1  “  pni2 


if  Y  .  <0,  then  p“ 
ni  nil 


Pni1  “  pni2 


and  for  6ni  <  0 


pni  1  “  -  Pnl2  *  Pni  +  l|Yni  +  I  6nil  >  /2  (23> 

pni1  ’  pni2  “  f^ni/2/pni1 

7 

The  range  for  the  various  types  of  roots  depends  on  the  magnitude  of  6ciaj[ 
for  each  of  the  values  of  6nj. 


For  S  >  0,  the  roots  are  real,  if  ip  ,  -  £  0  (meaning  Y  .  2  {  /), 

ni  ni2  ni  ni 

2  2 

i  a  ^  ^  A  i_i _  _  /  n  _  _  ^  \ 


thus,  if  6cla1  >  n  -  1 ,  then  snl3  *  s^  (Real)  , 


2  2 

and  if  6cla1  -  n  -  1,  then  s^  -  snll)(Real) 


The  four  roots  are  equal  if  6ni  -  0,  thus 


8cla2  -  -  4n  (n  - j  n2-l) 


( See  f i gs . 3 ) • 


The  various  solutions,  Fnik(r)»  depend  on  the  type  of  roots,  (real, 


equal  or  complex).  The  relation  between  6  ,af  and  P  is 

ci  i  cr 


70 


Complex  roots 


a  .a  D 

cl  1 


where  N  v  [see  Eq.  (6)]  Is  the  value  cf  N  per  unit  of  buckling  force, 
cl  cl 

and  D°  Is  the  relative  flexjral  rigidity. 

The  two  lines  appearing  on  Pig.  3  represent  points  of  singularity,  for 
a  given  n,  In  the  buckling  character  1st i c  function  [see  Figure  4],  The 
character  1st l c  function  Is  obtained  by  expanding  the  determinant  that 
yields  critical  loads.  This  determinant  results  from  requiring  a 
nontrivial  solution  to  the  set  of  homogenous  linear  algebrlc  equations  In 
the  undetermined  constants,  Akn|.  This  set  Is  obtained  by  using  the 
boundary  and  Joint  conditions. 


complex  roots  solution 
Real  roots  solution 
Two  eq-ial  roots  solution 
Four  equal  roots  solution 


P  ?  (11 

crs,  «  ( n  -  11  D./N  , 
i  l  ocl 

P  2  ( 1  ) 

crsl  «  -4n  (n-  n  -1 )  Dl/fT  , 
2  ocl 


The  characteristic  function  versus  the  critical  forces 
for  the  various  types  of  solution" 


% 


< 


c 


The  singularity  in  the  char ac ter  is t i c  function  is  due  to  the 
discontinuity  in  the  type  of  the  roots.  In  each  interval,  the  function  is 
continuous,  but  not  at  the  ends  of  the  interval  with  the  equal  roots.  At 
these  points,  the  function  is  approaching  zero  either  from  the  positive  or 
the  negative  side,  and  exactly  at  these  points,  the  function  can  assume  any 
value,  including  zero.  This  numerical  phenomenon  occurs  because  the 
solutions  for  complex  or  real  roots  converge  to  two  or  four  equal  solutions 
and  do  not  approach  the  appropriate  solutions  for  equal  roots.  The 
characteristic  function  is  calculated  through  the  determinant  of  the 
buckling  matrix.  The  cases  of  equal  solutions  means,  two  or  four  equal 
rows  in  that  matrix,  meaning  the  determinant  is  zero.  The  roots  for  n  «  0, 
1  and  n  £  2,  and  for  r  ®  £  £  «*  are  described  in  Table  1 . 

3.2  Types  of  solutions,  F i^ni(r) 

This  type  of  solutions  depends  on  the  type  of  roots.  In  general,  the 
power  series  solution  based  on  the  recurrence  formula,  Eq.  (17),  is 
correct,  for  all  four  Independent  solutions  (k  =*  1,2, 3, 4),  only  if  the 
roots  are  real,  i.e.,  they  are  different  and  the  difference  between  each  of 
them  is  a  non^-integer . 

The  solution  in  case  of  two  equal  roots,  (in  general)  equals: 


Wr)  '  ctFkni<r)  ln  r  Vi  r 


sknl*J 


J-0 


(26) 


where  P,(<ni  and  are  solutions  corresponding  to  two  equal  roots, 

C i  is  a  coefficient  that  can  be  obtained  through  the  recurrence 


procedure,  and  b^j  are  new  undetermined  coefficients. 


>3 


pltcltjr) 


V 


t 
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The  solutions  in  case  of  complex  roots  are  based  on  the  recurrence 
formula,  appearing  in  Eq.  (17).  The  two  solutions  are  determined  by 
considering  the  real  and  imaginary  parts  as  two  real  solutions  yielding: 


R 

s 


Fi1  - 

Snl1 r_  .  I 

r  [S  . ,  cos  (s  ,, 

nil  nil 

m  r)  h  snl2  stn  (sjj,  in 

Fi2  * 

R 

sni1[S  ,,  sin  (s1,, 
r  nil  nil 

in  r)  *  snl2  003  <3^,  in 

where 

R 

Sni1 

‘  0  ’  anl2  “l"2* 

, .M/2  M/2  R  N,  3ni1 

“1  w  >  r 

R 

Sni2 

”  ani2  air2  + 

,b,,N/2  N/2I  N.  3ni1 

“i  w  > r 

.  R  ,  I 

A  .  =  a  ..  +  la  ,, 

nij  nij  nij 

(27) 


(27.1) 


In  the  case  of  Bci  =  0,  the  root3  are  integer  and  the  second  solution, 
Eq.  (26)  is  equal  to,  Yn(r),  a  Bessel  function  of  the  second  kind  of  order 

n.  Using  the  power  series,  this  solution  becomes: 


F 


ni2 


C.F,  . 
1  Ini 


in  r  .  C2r"ln  r  .  ^  bnlJ 


j 

r 


+  s 


ni2 


(28) 


1 


I 


where  C1 ,  are  constants  calculated  through  the  recurrence  procedure, 

n  is  the  number  of  full  waves,  and  s  ,  „ 

ni2 

equations  (See  Table  1). 


t 


7  S' 


2-n  is  the  second  root  of  indical 


The  type  of  solutions  for  the  buckling  equation  depends  on  the  following 
parameters: 

(1)  Number  of  circumferential  full  waves:  n  -  0,1,2,...N 

(2)  The  value  of  ratio  of  the  non*-uniform  part  to  the  uniform  one, 

<  or  =  or  >  0. 

(3)  The  type  of  roots  of  the  indicial  equation,  e.g.,  real,  complex, 
equal  or  different  by  an  interger. 

The  various  solutions  for  the  parameters,  mentioned  above,  are  described  in 
Table  2.  Moreover,  the  various  solutions  are  shown  in  Appendix  B. 

The  solutions  derived,  using  the  power  series  method,  usually 
converge.  The  convergence  is  checked  either  by  using  tests  for  infinite 
series  [21]  or  by  substitution  of  the  solution  into  the  buckling  equations 
Eq.  [13).  The  criterion  for  convergence  using  the  second  approach  is: 

I  L  <FklnrJ)  I  S  ,0”6  <29) 

where:  L  is  a  differential  operator,  see  Eq.  (13), 

rj  is  the  radial  coordinate  at  the  boundaries  of  an  annular 
plate,  and 

Fkin(rj)  is  the  k?*th  solution  of  plate  i  different  from  zero. 

If  the  test  falls,  then  the  number  of  terms  in  this  series  must  be 
Increased.  The  results  obtained  by  using  these  tests  are  considered  to  be 
accurate,  when  compared  to  the  results  reported  in  the  open  literature 


(used  as  benchmarks). 
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2:  Solutions  for  fi 

i  root 

l  Solution 

srl**'i 


1*0  u2C <1*0J )2-n2: 

1  "  [(2+0,)2  -023C  <2+0,  >2+023 

(-1  >N/2UN  *rc<  J+0,-1  ‘  j-n2) 

....+N - ^ 1 r  > 

Tfi,(  j +u> t )  2+<i>2  3C  ( J+0J  )"-02  3 


2  2  2 
1-0  u  C<l-0  >  -n  3 
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S3=l+102 
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where: 
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V'fVz.  ‘  “  ‘ 

•1*' 
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0.,=  l+n  +13  +  C-<2n“+f3  */2>  (n  -1 )  3  = 
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_  D  2 
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n  J  n  n  Vz  Va 
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(  J4-1 )  2 — 0 2 )  }/  j  Dj  j  2 

2  2  2  '»  a  *»  ->  .-'i 

lDjl  =16  j  (j  -  v>  )  “+j  (  j"-4k>;-2^  >“■ 
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TABLfc  sS:  I 

~  "  •  '  1. 1 - 


n 

root 

1  Solution 

0 

Sl=2 

c  a 1 _  + 
F1  1  n2  *' 

N/2  N 

. - —  =  Jrt<u> 

N  -2  <-> 

T£<  j+2> 

*2“2 

t)F  (s) 
F2“1>. 

-  Y  .  (u  > 

5  =  5^  0 

s^=0 

•-> 

F_=Const . 

m 

F4=  In (r) 

Notesi 


JQ  <u)  -  Bessel  -function  o-f  -first  kind  and  o-f  order  0. 
Y^Cu)  -  Bessel  -function  o-f  second  kind  and  o-f  order  0. 
u=«T  r 

l 


Notesi 


».-2+n  F  =rn<l - - - - - - 

1  <2+n> A-n“ 


J  (u) 
n 


,  , ,N/2  N 
(-1 )  u 


2 

♦2+n ) -n  3 


=  Y  <u) 

s^s^  n 


J  (u)  -  Bessel  function  of  first  kind  and  of  order  n. 
n 


Y  <u)  -  Bessel  function  of  second  kind  and  of  order  n 


(Subscript  are  omitted  -for  simplicity) 


n  I  root  I 


Solution 


0^«sqrEl-abs(Rc^>at^  D 


Solutions  are  the  same  as  in  table  2,  but 


with  a  different  0^  . 


Solution 


same  as  in  table  2  with  n“0 

same  as  in  table  2  with  n=Q 

Same  as  in  Table  2  with  n=l. 

0.  =sqrCl-abs(R  .  )^2D 
1  Cl  1 

F,=r 

N  + 

F4=rln(r)^-j^0bjrJ 

The  series  coefficients  are  determined 
using  the  procedure  described  in  Appendix 
8.2.3.  1 


Fl=r 

F,*r  1  n  (r)+£^0 b^r  j  +  1 
F3»rCln(r>  32+ J^c  r J  +  1 
F4-rCln<r)33+J4d.rj  +  1 

The  series  coefficients  are  determined 
using  the  procedure  described  in  Appendix 
B.  2. 3.  1 
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The  series  coe-f  -f  1  ci  ents  are  determined 
using  the  procedure  described  in  Appendix 
B.2.3.  1 
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<0  (continued) 
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Notes (continued )  i 

r>  . 


Sol ut i on 


1  D®0 Ic‘4 j  (V*+20R>  +40^3 

o  o  «->  o  *  * 

abs(D. )  =R^+Ij0^ 

T’  =  l+n2+abs(fi  )V.2/2 

v  c  1  1 

T  =4n2+2n"1’abs  (G  .  )  'i‘+<abs(G  )  ?>-V4 

ci  x  Cl  1 

2  c"  V2 
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s^.=  l-0_-i  0. 

O  R  I 

1-0 

F^=r  CcosC0Tln(r) 3S. +si nC0T 1 n (r ) 3S_  3- 

o  1112 

s4=l-0R+i0I 

1-0 

F4=r  {cosC0Iln(r>  3S2-sinC0Iln  (r  ) 

-5 

The  varios  -functions  appearing  in  the 

solutions  are  the  same  as  above  but 

with  0R*-0R. 

-3 

TABLE  5 1  Solutions  for  fi  . 
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S3=2 


s4=0 


F^=Const. 
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For  n  -1<  B  . 
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F  1 =r 
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F  =r 


»,=  l+0. 


s4=l-0„ 


F_=r 


F4=r 


F3=rcosC02ln (r ) ] 


F^arsinCti^ln  (r  >  ] 


For  n 


S « - 1 +0 


Notesi 


F4»rln (r ) 


0  =<l+n2+ft  .  /2>+<4n2+2n2R  .+Q2./4)Z 

1  ci  ci  ci 

0_=<l+n2+ft  . /2>-<4n2+2n2ft  . +fi2. /4>^ 

2  ci  ci  ci 


The  solutions  are  the  same  as  above  but  with 

I  ft  .  |  instead  o-f  ft 
i  ci l  ci 


.  NUMERICAL  ANALYSIS  AND  RESULTS 


( 


C 


A  computer  program  was  developed  for  the  analysis  of  multinannular 
circular  plates,  for  various  boundary  conditions  and  types  of  loading.  The 
program  uses  double  precision  and  the  maximum  number  of  terms  in  each 
series  is  eighty.  The  "differential  equation"  convergence  condition  [(Eq. 
29)]  is  checked  for  each  independent  solution,  F^itr),  each  coordinate 
of  the  adjoining  boundaries,  r j .  The  various  solutions  are  obtained 
through  the  following  four  subroutines. 

a)  Solution  for  uniform  stress  field,  meaning  Bci  «  0. 

b)  Solution  for  real  roots,  but  distinctly  different. 

c)  Solution  for  two  or  more  equal  roots. 

d)  Solution  for  complex  roots. 

A  relation  exists  between  the  various  subroutines.  The  solutions 
obtained  using  subroutine  (b),  for  integer  roots,  converges  to  the  one 
determined  by  using  subroutine  (a).  The  solution  for  nearly  equal  roots, 
obtained  by  subroutine  (b)  for  real  roots,  or  subroutine  (d)  for  complex 
roots  with  equal  real  parts  and  very  small  imaginary  parts,  do  not  converge 
to  the  solutions  obtained  by  subroutine  (c). 

Numerical  results  are  generated  and  discussed  f or  a  two  part  plate, 
annular  and  circular  (inner),  with  constant  thickness  only.  The  main 
parameter  influencing  the  buckling  force  are: 

1)  The  ratio  of  the  radii  ri/r0,  inner  radius  to  outer  one. 

2)  The  ratio  of  the  moduli  of  elasticity,  E2/E1 ,  inner  part. 

3)  The  number  of  waves  in  the  clrcumf erential  direction  (buckling  mode). 

*0  The  location  of  the  loading,  applied  at  the  outer  radius  as  opposed  to 
the  inner  one. 

5)  The  boundary  conditions. 
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In  order  to  get  an  understanding  of  the  influence  of  the  various 
parameters,  a  few  cases  are  analyzed,  and  some  simple  parametric  studies 
are  performed. 

4.1;  Example  No.  1  (Simply^-supported) 

The  first  example  is  a  two-part  plate  loaded  and  supported  at  the 
outer  boundary,  see  Fig.  5.  Critical  (buckling)  loads  are  calculated  for 
the  following  parameters: 

v-j  -  \>2  =  1/3  (n«  0,1  modes) 

t-j  =*  t^  —  1.0 

E2/E1  =  0  (Hole),  0.01  (cir.  part  very  flexible), 

0.1,  0.5,  1  ('uniform), 

2,  10,  100  (cir.  part  very  stiff) 

r-\/r0  -  0.0,  0.1  .  0.9  (Dec.  =*  0.1) 

The  analysis  using  the  above  data  covers  the  whole  spectrum  of 
rigidity  ratios  starting  with  a  single  annular  plate  (E2/E1  *  0),  through  a 
plate  with  a  very  flexible  central  part  E2/E1  =  0.01  to  a  plate  with  very 
rigid  inclusion  (E2/E1  -  100). 

Curves  of  critical  loads  versus  the  ratio  of  the  radii,  r^/rg,  for 
different  (E2/E1)  values  for  n  =  0  (axisymmetric)  and  n  =  1,  with  v  »  1  / 3 
are  presented  on  Fig.  5  and  6,  respectively.  From  these  two  figures  the 
following  conclusions  can  be  drawn: 

1)  The  lowest  critical  load  corresponds  to  the  first  mode,  even  for  large 
values  of  E2/E1 . 
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Fig.  5.  The  buckling  Force  versus  the  radii  ratio  (ri/r0)  For  various 
rigidity  ration  with  n  -  0  (axisymmetric  mode)  and  a 
s i mpl y - s uppor t ed  plate.  The  plate  Is  loaded  at  the  outer 
boundary. 
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2}  The  critical  load  decreases  with  increasing  value  of  r^/r,-,  for  E2/E1  < 
1.0.  On  the  contrary,  for  E2/E1  >  1.0  the  critical  load  increases 
with  increasing  value  for  ri /rg. 

3)  For  stiff  inclusions  (E2/E1  >  1),  the  buckling  mode  is  such  that  both 
parts  of  the  composite  plate  participate  (the  inner  complete  circular 
plate  and  the  annular  part),  and  buckling  seems  to  be  triggered  by  the 
annular  part. 

4)  For  flexible  inclusions  (E2/E^  <  1),  the  buckling  mode  depends  on  the 
values  of  both  ratios,  E2/E1  and  r^/ro.  For  large  values  of  E2/E1 
(definitely  for  0.5  £  E2/E^  <  1),  the  mode  behavior  is  similar  to  that 
corresponding  to  the  case  of  E2/E1  >  1  (See  3)  above).  On  nhe 
contrary,  for  extremely  small  values  of  E2/E1  (it  is  checked  for 
values  of  0.1  and  0.01  and  the  observation  probably  is  true  for  some 
E2/Ei  values  in  the  range  of  0.1  <  E2/E1  <  0.5)  the  mode  behavior 
depends  on  the  ratio  of  the  radii  (ri/rg).  For  E2/E1  =  0.1,  as  long  as 
r-|/ro  is  smaller  than  approximately  0.55,  the  buckling  mode  seems  to  be 
triggered  primarily  by  the  annular  part  and  the  critical  load  is  a 
little  larger  than  the  one  corresponding  to  E2/E1  =  0  (pure  annular 
plate).  On  the  other  hand  for  ri/rQ  *  values  larger  than  0.55,  the 
critical  load  is  smaller  than  that  of  the  pure  annular  plate  and  the 
buckling  mode  is  triggered  by  the  inner  and  more  flexible  part.  It 
appears  that  for  E2/E1  =*  0.1,  the  corresponding  curve  approaches  the 
value  of  0.4282  as  ri/r0  approaches  unity  (a  complete  circular  plate 
with  bending  rigidity  equal  to  one  tenth  of  that  of  the  case  E2  -  Ei). 

5)  The  graph  for  a  single  annular  plate  (E2/E-{  »  0),  is  similar  to  the  one 
appearing  in  Meissner's  paper  [10].  The  only  difference  is  in  the 
result  for  the  homogeneous  plate  (E2/E1  »1).  The  present  result  for  the 
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critical  load  la  4.282D /r^  and  not  4.20D/r2  as  predicted  by  Meissmer. 
The  present  result  corresponds  to  v  *  0.333  and  Meissmer’ 3  to  v  -  0.3. 
Thus,  the  difference  is  a  measure  of  the  Poisson's  ratio  effect. 

It  appears  that  the  influence  of  the  Poisson's  ratio  is  significant 
for  annular  plates  with  either  small  or  large  values  of  the  ratio  for  radii. 
For  small  values  (r-]/r0  <  0.4),  the  higher  the  Poisson's  ratio  the  stronger 
the  configuration  (against  buckling).  For  r^/r0  -  0  (full  circular  plate), 
the  difference  between  the  critical  loads,  corresponding  to  v  -  0.333  and  v 
»  0,  is  of  the  order  of  over  25?.  This  difference  diminishes  (for  the 
annular  plate)  as  r^/r q  approaches  the  value  of  0.5.  For  higher  values  of 
ri/rQ,  the  influence  reverses  direction  and  it  increases  with  increasing 
ratio  values.  Note  that  for  r-j/rQ  -  0.9,  the  critical  load  for  \i  *  0  is 
higher  than  that  for  v  «  0.333  by  approximately  10?. 

4.2  Example  No.  2  (clamped) 

In  the  second  example,  the  influence  of  the  boundary  conditions  is 
checked,  through  the  analysis  of  a  clamped  plate. 

In  this  example,  the  clamped  plate  is  loaded  at  the  external  boundary. 
The  critical  load  curves  for  a  single  annular  plate  are  presented  on  Fig.  8 
and  for  the  two-part  plate  on  Fig.  9.  The  results  are  determined  for  the 
following  parameters: 
vi  ■  \>2  *  1^3 
C1  ■  fc2  -  1 .0 

n  -  0,1, 2, 3.4  (wave  number) 

e2/V-\  -  0,0.1,0.5,1,2,10,100  (n  -  0) 

ri/r0  -  0,0.1  ,  0.8 
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The  critical  load  curves  versus  the  ratio  of  the  radii,  for  a  single 
annular  and  for  wave  numbers,  zero  to  four,  are  presented  on  Fig.  8.  The 
curve  for  the  axisymmetric  mode  is  exactly  the  same  as  the  one  appearing  in 
Meissmer’s  paper  [10],  The  results  appearing  in  Ref.  [10]  are  limited  to 
r/rQ  £  0.563.  for  ratios  greater  than  0.52  the  higher  buckling  modes 
govern  buckling  which  is  in  contrast  to  the  previous  example  (simply 
supported  case).  On  Fig.  9,  the  critical  load  curves,  corresponding  to 
different  E2/Ei*values,  n  =  0  (anisymmetr ic ) ,  and  v  -  1/3  are  presented. 
In  this  case,  Fig.  9,  the  following  conclusions  can  be  drawn. 

1)  For  E2/E|  >  1  the  critical  load  increases  with  increasing  r^/rQ  - 
values,  while  for  E2/E1  <  1  the  picture  is  a  little  more  complex. 
First  of  all,  spot  checking  of  these  results  (E2/E1  >  1)  showed  that 
the  critical  condition  is  governed  by  axisymmetric  behavior  (n  =  0). 
It  is  also  observed  that,  as  the  middle  section  (the  complete  circle) 
becomes  stiffer  the  critical  load  increases  more  rapidly  with 
increasing  r^/rg  '-values. 

2)  For  E2/E1  <  1 ,  we  observe  several  things.  First,  as  long  as  the 
stiffness  of  the  middle  section  is  only  two  to  three  times  or  less 
smaller  than  the  stiffness  of  the  annular  (outer)  section  the  response 
is  similar  to  that  of  the  simply  supported  case,  and  there  is  a 
decrease  in  value  of  the  critical  load  as  the  ri/rQ*value  increases. 
In  this  case  the  buckling  mode  is  such  that  both  parts  (inner  and 
outer)  participate  equally. 

On  the  other  hand,  when  E2/Ei  is  zero,  we  have  already  seen  (Fig. 
8)  that  the  critical  load  increases  as  rj/r 0  Increases,  and  the  mode 
becomes  non^axisymmetric.  For  E2/E1  very  small,  or  for  inner  section 
stiffnesses  several  times  smaller  than  the  outer  section  stiffness,  it 


Is  observed  that  for  low  values  of  the  ri/rg  ratio  the  behavior  is 
similar  to  that  of  E2/E1  -  0.  This  means  that  there  is  initially  a 
decrease  of  the  critical  load  as  r^/r0  increases;  a  minimum  value  is 
reached;  and  then  there  is  an  increase  (see  curve  corresponding  to 
E2/E1  -  0.1).  As  the  value  of  ri/rg  further  increases,  there  is  a 
distinct  decrease  in  the  critical  load.  The  explanation  here  is  that 
in  the  former  case  buckling  is  governed  primarily  by  the  annular 
(outer)  section  while  for  the  latter  case  buckling  is  governed  by  the 
inner  and  less  stiff  part.  Clearly  in  both  cases  (for  the  full  range 
of  r-|/r 0^ values) ,  the  presence  of  a  very  flexible  inner  part  weakens 
the  annular  plate.  Note  that  the  curve  corresponding  to  E2/E1  *  0  is 
higher  than  that  corresponding  to  E2/E1  -  0.1. 

4.3  Example  No.  3  (Load  at  r  ■  n  in  the  Inward  direction) 

In  this  case  a  simply**supported  plate  loaded  at  the  inner  radius  is 
considered.  The  uniform  radial  load  is  in  the  outward  direction.  The 
curves  for  the  following  data  are  presented  on  Fig.  10. 
vi  •  V2  *  1/3 

ti  <■  t2  *  1.0 

E2/E}  -  0  (Hole),  0.05,  0.1,  0.15,  0.2,  0.3 
rj/r0  -  0.2,  *••-»»,  0.8  and 
n  -  0  (Axisyum/netric  mode) 

The  curves  on  Fig.  10  correspond  to  a  single  annular  plate  with  a 
flexible  inclusion.  The  critical  load  for  a  single  annular  is  very  large 
for  ri/rg  <  0.2,  since  the  inplane  stresses  are  very  small  in  most  part  of 
the  plate  and  are  very  large  near  the  point  of  application.  For  ri^ro  > 
0.2  the  stresses  are  distributed  more  uniformly  within  the  entire  plate. 


Fig.  ii:  The  total  buckling  fore*  (Pcrrl>  ',«rsja  ri/ro  for  various 
rigidities,  with  n  -  0  and  slapljr  supported  plats,  loaded  at  the 
inner  Joint. 
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In  the  case  of  the  plate  with  the  flexible  inclusion  the  buckling 

\ 

force  becomes  small  around  ri/rg  p  0.55.  The  critical  force  increase,  for 
all  values  of  E2/E1  ,  as  the  radii  ratio  increase,  since  the  tension  in  the 
central  part  is  governing.  The  force  for  all  ri/rg  increases  with 

> 

increasing  of  the  rigidity  ratios. 

The  "total"  force,  pQr  r^ ,  versus  radii  ratio  ri/rg  is  described  on 
Fig.  11.  In  general,  the  trends  of  the  "total"  force  appearing  in  the 

►  figure  and  the  critical  force  described  in  the  previous  one  are  the  same. 

The  "total"  buckling  forces  for  ri/rg  £  0.2  are  bounded  and  don't  go  to 

infinity  as  can  be  deducted  from  Fig.  10. 

5.  CONCLUSIONS 

A  buckling  analysis  of  a  multi-annular  plate  is  described.  The 
stiffness  method  Is  used  to  determine  the  primary  state.  The  buckling 
*  equation  is  solved  rigorously  using  the  power  series  method.  Convergence 

is  proved  either  by  the  infinite  series  convergence  test  or  by  the  accurate 
satisfaction  of  the  differential  equation  condition.  The  versatility  of 
^  the  proposed  solution  is  proved  by  the  solved  examples.  The  numerical 

examples  reveal  that  the  influence  of  the  different  annular  parts  and  of 
the  central  part  may  or  may  not  improve  resistance  to  buckling.  The 
Influence  of  the  various  parameters,  ri/rg,  E2/E1 ,  v,  n,  is  discussed  to 
3ome  extent  for  a  two-part  composite  plate.  A  complete  parametric  study  is 
recommended  to  understand  the  influence  of  the  various  parameters  on  the 
buckling  behavior  of  multi^-annular  plates. 

The  proposed  solution  is  also  applicable  to  the  buckling  analysis  of 
thin  plates  when  resting  on  an  elastic  foundation,  and  even  to  the 


vibration  analysis  of  multi*annular  plates 
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Appendix  A:  Derivation  of  the  Buckling  Equation  and  Inrplane  Stiffness 


Coefficients. 

A. 1  Derivation  of  the  equilibrium  equations 


The  nonlinear  kinematic  relations  in  polar  coordinates,  are  [22] 
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where:  ey  f~  strains  of  the  midplane  (i,j  -  r,0) 

u,v,w  ri  displacement  components  of  midplane  points  in  the  x,y,  and  z 

directions,  respectively  and 
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differentiation  with  respect  to  i: 


The  equilibrium  equations  for  isotropic  material  (see  Fig.  A.1)  without 
body  forces  are: 
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f  flexural  rigidity  of  the  plate. 


Nrr’  N00 *  Nr0  "  Stress  resultants  (for  sign  convention  see  Fig.  A 
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DcforMtlont  of  a  differential  Plate  Area 


< 
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The  Laplacian  operator  V2,  operating  twice  on  w,  yields 


4  2  1  1 

V  w  -  w,  +  —  w,  +  — -  (*  w,  +  2w,  ..  )  +  —  (w,  u  2w.  .  ) 

rrrr  r  'rrr  2  rr  ’  00rr  3  *r  90r 

r  r 


4  (4w,00  +  w’  9009 ^ 

P 


(A. 3) 


Next,  the  hackling  equations  are  derived  using  the  perturbation  approach. 

The  stress  resultants,  and  transverse  displacement,  in  the  buckled 
state,  are  given  by: 


rr 


00 


o  * 

N  +  N 
rr  rr 

o  * 

+  Nqq 
00  00 


(A. 4) 


r0 


o  * 

AT  *  N  Q 
r0  r0 


o  * 

W  -  W  +  W 


where:  (  )  s*  primary  state  parameters 

# 

(  )  **  additional  small  parameters  need  to  reach  a  buckled  state. 

Substituion  of  the  quantities  of  Eqs.  (A. 4)  into  Eqs.  (A. 2)  and  some 

algebraic  manipulations  yield  (note  that,  j  j  <  <  |  N  I  ancl  w°= 


(r  N*  )  +  N*  *  N*  «  0 

rr  ,r  r0,9  90 


N*  +  (r  N*  ) ,  +  N*  =•  0 

90  0  r9  r  r0 


(A. 5) 


4*1o*  o*  o*  o* 

DV  w  -  -  r  N  w,  ),  +  (N°  w,  Q/r ) ,  .  +  (N°  w,Q),Q  +  (N  °w,  J, 

r  rr  'r  r  90  0  0  r0  0  0  r0  9  r 
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The  expressions  for  the  bending  moments  and  transverse  shear  forces  are 
[19]: 


M  »  '-D  [w,  +  v  (  —  w,__  +  ■— w |  )] 

rr  ’rr  2  90  r  *r 

r 


M.„  =  -D  [w,  /r  +  — -  w,„„  +  vw,  ] 
00  ’r  2  00  *rr 

r 


Ma”*-(1*v)D[—  w,  . - -  w ,  a  ] 

r0  r  ’r@  2  0 

r 


V  -  —  [ (M  r),  +  M  .  .  M  ] 

rz  r  rr  r  r0,6  80 


-  D  (V  w). 


(A. 6) 


VQ  -  -  [M,aa  +  M  „  r  +  2M„]  =  ^-D  -(V2w),„ 
9z  r  ’00  r0,r  r0  r  ’8 

and  the  edge  effective  shear  forces  are: 


Q-V  +  -  M  „  Q 
r  rz  r  r0,0 

Q  =  V  +  M 
0  0z  r0,r 

A. 2:  Primary  state  analysis 

The  stresses  and  the  deformations  at  the  adjacent  joints  of  the 
annular  parts,  are  calculated  by  using  the  stiffness  approach.  The 
formulae  for  stresses  are  taken  from  Ref.  [20]. 

A. 2,1  Stiffness  Coefficient  [3&j] 

The  stiffness  coefficients  are  determined  separately  for  the  circular 
and  annular  parts. 
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The  relation  between  the  off-diagonal  terms  is: 


ri  3i-1  i  “  riH  3i  i-1 


and  the  stressses  and  deformation  in  the  plate  are: 


(r  £  r  £  r  ) 
v  i  i'-l 


N°  (r) 
rr 


[(1  +  )  + 


8i(i-Vri 


u 


r^Bj  *•  1)(1-v‘) 


2  C(1  +  v1)(-ajl+  ) 


8i(1~vi)ri 


("8iV  'J: 


2  2 
r .  r. 


u(r)  -  [u.  B,  (1  -  (  —  )  +  u  (  *l)  ] 


Appendix  B:  Various  Solutions  of  the  Buckling  Equation 


B. 1  Solution  of  the  buckling  equation* 


The  buckling  equation,  after  using  separation  of  variables  reads: 


rVIV  +  2r3Fm*  r2(1  +  2n2+  e\a2  <*■  r2a2)  F11  +  r(1+2n2+  8  a2  +  r2a2)  F1 

Cl  1  1  Cl  1  1 


-F  [n2(i»-n2  +  8  a2  +  r2o2  ) 
Cl  1  1 


(B.  1  ) 


where 


1  D, 


and  8  .  * 


cl 


cl  N 


We  first  assume  a  power  series  type  of  solution 


F  (r)  -  I  A  r 
j-0  J 


j+s 


(B.2) 


Substitution  of  this  solution  Into  Eq.  (B.1)  (after  some  algebrlc 
manipulations)  yields 


I  A  rj+3{(j  +  s  -  I)4-  (2  +  2n2*  S  .a2)  (j  +  s  H )‘ 
j-0  J 


,2„  2- 


[(n-1)2  +  8cia2  (1«n2)]  }  +  \  A  rj+S+2a2  [(j  ♦  s)"“  n6]  -  0 


(B.3) 


j-0 


*For  simplicity  the  subscripts,  nl  or  1  or  k  are  omitted,  throughout  this 
appendix. 


The  Indiclal  equation  reads  (j*0): 


(sM ) (2  +  2n2  +  Bcia^)  (s*-1)2  +  [(m-1)2+  8cia2  (l^n2)] 


and  the  roots  of  the  quart .tc  equation  are: 


02  H 


)?  [ (1+n2  +  —  {  —  )1  +  \L2  +  2nV  ,<*2  +  --V-  ■  * 


(sM)^-  L(1+n‘"  ♦  )J  -  ^n2  ♦  2n28cla2  ♦ 


T-  4  *  ,  (B.5) 


S1 ,2  *  1  1 


\|y~+~6^ 


93,4  "  1 


(B.6) 


where  Y  Is  equal  to  the  terms  appearing  In  first  brackets  and  6  equals  to 
the  terms  appearing  under  the  square  root  sign  In  Eq.  (B.5).  The 
recurrence  formulae  Is: 


af[(j+3"2)2*-  n2] 


a?[(j+s'-2)2  *-  n2]  A 


B.2 


Solutions  of  the  Differential  Equation 


B.2.  l:  >P 

The  solutions  for  a  positive  (3^  are  derived  next. 
B.2. 1.1.  n=0 


For  this  case  the  roots  ares 


1  -  (  i  4  ^  e//')  ^  -  \  ~  • 

*2>'  ^ 


S  4*0 


where 


f ^  ••  f*  .  '•  ^  s-** A  • 


The  recurrence  formula  is: 


^  ( i-*  s>-2^;  N;-2.  _ m  (j* 

1  '] 


T 


4^*1  *0 


('fc»AO>\ 


The  general  solution  is 


P  k  sf,_dN^ _ , 

F=  t  . 

ion  for  the  various  roots  are:  v)  ^ 


The  solution 


■**) 


-■  r  \  k  '-^»v  |  1  V<^ _  T  -  •  * - ‘ 

>-•*>;) ^ 


10S 


«4 

Vs 


^  =  r 


oc ;  [  C  j+s^-  »!)  -  n  "\  \.(s-a)  ^  (>  \)  b». ;-; 

3 

where  D.E<Fv<)=  Eq.  (s»l^with  F=Fy^.  Thus^  the  result  is:  /" 

^.E.  fc)  >  Ck-0>s-o V3  *-£ri j 

Equating  the  coef  f  i  c  i  ents  of  the  various  powers  of  r  to  zero  yields: 


r  '■  ^‘*'3  =4. 

->  \  Z*;* 


W,.Z*£ _ 

li^^Vyrj^y- 


‘"(rjv'  uh  1  ^  Urv*)V]  -< 

l  ■  -«?Ui»s».-x^rrt  ■_  1 


hence  the  fourth  solution  is: 


^  ^  L'hOr)  +  '  v  ^ v  ** 

r2  ] 


B. 2. 1.3:  n  >=2 


When  n>“2  the  roots  are: 


< 


where: 

f  * 

,  x. 

(2.o\  AoMoV,') 

(■&.«. f) 

"V;* 

vw 

r 

The  recurrence  formula  is: 

vJote  that  the  detailed  solution  is  described  in  article  8. 


B. 3. 1.4:  A  qeneral  solution  for  a  case  with  complex  roots 
The  complex  roots  are  described  by! 

S»(X,k)  =  ^  \SI(XNV)  (  ^ 

and  the  recurrence  formula  in  general  is: 


*  o*v\  j 


where: 


After  some  algebraic  manipulations  Me  obtain: 


where: 


•*  v  ^X-^CnrCL 


T\RA= 

T\^2  «.  ^vr^_-t\1^ 


^U-2.-  2'T\Ki  T\R£_ 


The  nominator  is: 


£Vrvow\  •  i„  \  *2> -•£*)*■—  T")'  s  pV rsorv> V 


where: 


JXrNvVV^  d 


TT\ 


-  (^s^-«aT--2>I 


1  »  ^(rtS^fs.l') 


The  general  recurrence  formula,  using  complex  notation,  becomes: 


VvV;M^r 


where: 


^\C.  ~  ~  - - 

^  *fc>  ^  wwe-wiT.”1" 

1  K^e-wvT  •<0«aa»v«.R  — 

&**“ - — - ; - ^ - 


*29>.\) 


For  5x85!,-,  ,it  can  be  easily  shown  that  the  recurrence  equation  is: 


Hence  the  general  solution  is: 


1=.  tr 

V^>oV  <£- 


t'  -  '  ^■at)  d  \  r  * - ♦  (-0  (i(  r- )  If  (&■*•» 

♦);  a.  J 

[  K  -  ^  °*  ;**•*•-  -  .  .-*(-0  (<*;/)  TT  (Vfe-v 

U  V  /  |-a-  1 


**  ^t-Px+F-j  and  F1=Fl-F3  with  A<~.=B,,= 1  /  2  then,  after  some  algebraic 
manipulations  we  get 


V’  (si  Uv(»l)  —  ^3,  *sAr»  f  SvV^(r->')J 

^J=Y"  Csi Vr*  |V\) ^*C 


S>'=  4  *“  -  *  -  -  +(-\?^Cjitr f 


where: 


w 


8.2.2:  3C,  =0 


The  various  solutions  tor  Gc, =0  are  determined  next. 

B.2.2. 1  n=0 

The  roots  are:  S,3S4-2. 

-  SA -  O 


Ml 


th Y=1  ,C0.  The  recurrence  formula  is  (for  sf=U)  t 

Pv -X  =  _ N  . 


Hence  the  general  solution  is: 


V 


I  - 


W.rQ 

yo  J 


3*^ 


V&\ 


^2>-  A. 

F^s  V^OO 

where  J  (  ,r>  and  Y  <  ,r)  are  Bessel  functions  of  the  first  and 
second  kind,  respectively  of  sero  order. 


< 


B.2.2. 2:  n=0 


The  roots  are: 


It  can  be  shown,  very  easily,  that  the  various  solutions  are: 


where  J nWir)  and  Yr,(i^1r)  are  Bessel  -functions  o-f  the  first  and 
second  kind  respectively  of  order  n. 


The  solutions  -for  n=0  and  n >=  1  are  determined  in  the  next  chapters. 
B. 2. 3. 1 :  n=0 


The  roots  are:  »  2 


where 


N2- 

' V 


For  abs  (Bc- 1  )&t,.  =  1  the  roots  are: 

o 


■Afl 


The  first  two  solutions  appear  in  art.  B.3.1.1.  The  last  two  are 
derived  next. 

The  recurrence  formula  is: 


77^  4  , . 


1  V 

Hence  for  s-.  =  l  we  have; 

L1  r  -  ^  n 

i  r 

The  third  solution  then  is: 


_  (i 

1 

r  . 

>  '  _  S 

C^- 


>  • 


-  - - ^13_J - 


The  fourth  solution  is  assumed  to  have  the  following  form; 

(X> 

\’3' 


C^U(r^  -v  22.  S'r  ^ 


f’&A'i) 


Substitution  of  the  solution  into  the  D.E  Eq.  (B.l),  yields: 

-  C,^  < ^5.  yCyS'-i.') 

'K  (V’Vh-,)  _2.^*  fc,! 


(iW^')  =vZ_y'  ^  w  -V 

^eTX  ^  ^  ^  ^  ^  ^  ^  T ^ 
fV'nQVyv  t 


where: 


Equating  the  coefficients  of  the  various  powers  of  r  to  zero 
yi  el ds: 

O  O  -J 

C\ic^  •  v)  Cl«-C  t  -V-  <=  -<=* 


we  choose  : 


C:*Vq  'f  •  < 


V^.^irT^  :  f  SS4\^  "Sv,  C2.  -t  ^  2>^ 


4-  ts^  ^<-« 


- ■  (s=>.4<) 


where  DAI  equals  the  terms  in  the  large  brackets  of  Eq.  (8.44) 
The  recurrence  formula  for  b is: 


L  •  _  -tA*\  -^=>(-2.  ^^OTfN 


B. 2.o.2:  n=l 


The  roots  are: 


SA=  =\4^>- 

^  i- (A'  U Vv 


\ 


where 


The  third  and  fourth  solutions  are  the  same  as  those  derived  in  art 
B.2. 1.2. 

The  first  two  roots  depends  on  0t  .  If  abs(B,:1)  ,<4  then  those  two 

roots  are  real  and  the  solutions  are  the  same  as  those  derived  for 
positive  mr;1  ,  (see  art.  B.2.1.2). 

For  abs(Brf )  ,=4  the  first  two  roots  are  equal  to  1,  namely  all  four 

roots  are  equal . The  four  solutions  are: 


e  trlxv(r)  ^ 

1~  1 


The  various  coefficients  bj ,  c ,  ,  dj ,  are  determined  by  the  same 
procedure  appearing  in  the  previous  article. 

For  abs  (Be- 1  )o<  <  !:  4 ,  the  first  two  roots  are  complex  and  equal  to: 


The  general  solution  for  this  type  of  roots  appears  in  art. 
B.2. 1.4. 


B.2.2.3:  n  >-2 


The  roots  are:  S> 


S(«  \  1-^ 

v  '  '  -^v 

&.»  \  *  \  y1'-  \ Vj-  ?  \  -v-v. 


For  abs  (I3c=t  )°<  i  <2.'S  (see  Table  l>;<i  >0  f  is  a  real 

complex.  Hence  the  roots  are: 


number  and' 


(^o.  5o] 


where: 


r  r'* 

S''5- 

^ 1  '  ■*  'Ki  - '  ^v,\  c*N  / >~ 

{“  =  f  2.^-  ^  C/-A  Y  -  A’-  0 


The  solutions  for  these  type  of  roots  are  given  in  art.  b.^.l.j- 
That  solution  also  covers  the  range  0<  abs(R(;:;1)  x  <2  . 

If  abs(Bcl  then  is  negative  meaning  the  roots  are: 


^.SV\ 


where: 


S»2=  V-V^Vv*  'K  *  \  -  v  AY 
V^.-  :  y» 

^  ^  2.  :  t  ^  w^r 

'■  o rc\«  b"  =  ^ 


(&&v0 


\  *C  1  2-  / 

’  ortV«^-=  ^  ,/yH 


and  5l=57  and  s^s*  (The  (  )  “■  means  ct  complex  conjugate  number) 
The  general  solution  -for  this  type  of  roots  appears  in  art.  B.2. 


In  this  case  the  uni-form  part  o-f  the  stress  -field  is  zero  and  the 
differential  equation  (Eq.  B.l),  becomes  an  Euler  type.  The  solution 
for  this  type  of  equation  is  given  by  functions  then  by  power 
series. 

The  indicial  equation  is  the  same  as  for  the  power  series  solution, 
Eq.  <B.4>.  The  various  roots  are  described  in  Table  1.  The  solutions 
for  the  various  types  of  roots  is  derived  next. 

The  roots  in  general  are: 


z  A  . 


B.  2. 4.  1 :  n=0 ,  >p 

The  D.E.  is  equal  to: 


and  the  roots  are:  <>,<=  ?“ 


\\\ 


The  general  solution  is: 


F.C/V'*C,^ 


cvv  +  cs 


B.2.4.2:  n»0,  (3c1<q 


The  roots  are: 


^a-2.*;  ^c° 


Hence  the  solution  is: 

*.  ->  Cjjr\  Ch  Mf) 

In  terms  of  tri gonometri c  functions,  the  solution  becomes 


C^Y*  r^*f  C \  -*•  C^V”  +  Q_s^  f 


B.  2.  4.  3:  n>0,  I3el  >0 
The  roots  in  general  are: 


-  \  +  ^  * 

\-{  ^V'*- 
**■*,-'  \*\ 

tc 

ll 

\ -  A^-V7*- 

where: 

V- 

^  =  A *  /(A 


The  roots  s3  and  s.*  are  real  if: 

r*>T 

Expressing  this  in  terms  of  n  and  I3C,  yields: 

Thus,  for  \r?— a  >  the  solution  is: 


and  for  Y\*- 1  <fe»e.',  we  get: 


Mc,^C,r^  C,r<F\  ^  ]  M 


When  we  get: 

>cizr^TJ,  ^  aHU(r)1 


■fc>.0 9*1 


\  22. 


B.  2. 4. 4:  p. >0 ,  ra.-t<0 
In  this  case  the  roots  ares  \  * *\  '*• 


‘Tv  ' 


,r'-  * 


and  Ip  and  C  are  the  same  as  in  article  B.2.4.3,  but  with  a  negative 
,  The  roots  can  be  real  or  complex  quantities.  The  roots  st  and 
s are  real  i  -f  s 


The  roots  s,  and  s«  are  real  if: 

o 


After  some  alqebraic  manipulation  and  substitution  o-f  n  and  fiM,  we 
get  j 


which  is  also 
The  solution 


l-Y^>  1  1 

r-\ 

also  correct  -for  IT*-'  > 

ion  -for  \  -  T'v  >  \  ^  IS! 


c  f  nTJ71  -{77^ 

F=r[t/T’  -<v'~ 

and  -for  the  solution  becomes: 

f.  e.3WnU)j4Vti\«NM]-t  Ci 

-c„^  1  M 

If  Y<*  the  solution  equals: 

T-^c.r'T^l  Cl-‘rW,*^C5^(r^fW.,l 


Appendix  C:  The  characteristic  equation  for  a  circular  plate  composed  of 

annulars  parts 


The  characteristic  equation  is  obtained  by  employing  the  boundary 
conditions,  and  compatibility  and  equilibrium  at  the  inner  joints. 

The  conditions  which  must  be  satisfied  at  the  inner  joints  are: 


at  r  «  r, 


W  »  w 

11  12 


3u,  ,  3w,  2 


3r 
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3r 
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(C.1  ) 


rrl  1  rr1  2 


Q  -  *  Q 
rr  rr 

11  12 


where 


„  .  [  if*  .  i  (|s  ,  iifs  ,j  D 

rr  „r2  r  3r  r  3e2 


(C.1 .1 ) 


_  r  3  n2  1-v  3  r  3  ,  1  3w  >,i  n 

Qrr*  ‘  3r  V  w  r  38  ^  3r  (  r  30  )] ^  D 


The  condition  at  the  boundaries  can  be  any  one  of  the  following:  (See  [9]) 


3w 


(1)  Clamped  edge:  w  -  0  ,  -  0 


(C.2) 


(2)  Simply  supported  edge:  w 


0 


(3)  Free  edge  subjected  to  uniform  compressive  force  p0: 
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32w 


+  _v  j-  _3w  +  _l_  d  w  I 
?  r  1  3r  r  ‘ 
3r  36 


i  >>2 

1  3  w 


rr 


(C.  4) 
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(V2w)  +  1_  r  i_  rl  3w)]+_°iw 

K  w;  r  30  1  3r  V  36  JJ  D  3 r  * 


(4)  Free  to  deflect  but  not  to  rotate: 


3w 

3r 
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(C.5) 
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rr 


J_  (V2w)  + 1-  r  i_  r  1  }] 

3r  K  ’  r  30  1  3r  1  r  36  11 


Po  3w 
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A  support,  at  an  inner  joints,  is  characterized  by 


r  -  rk: 


wk  kM 


k  k 
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3w 


k  k"1 
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(C.6) 
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In  order  to  derive  the  characteristic  equation,  we  start  with  the 
inner  part,  which  can  be  either  a  hole  or  a  continuous  circular  plate.  The 
case  of  a  plate  with  a  hole  is  described  by  boundary  conditions  at  the 
inner  edge  of  the  annular  plate. 


k 


» 


( 


The  case  with  an  Inner  circular 
plate  will  be  discussed  next.  The 
deflection  equation  for  the  circular 
part  is  described  using  SCi  -  0 
(See  Appendix  B)  for  various 
n*-  values. 


v> 


C,„VU>  *  C2nVu>  *  V 


C,  un 
4N 


(C.7) 


where  u  »  a^r. 


For  F(r)  to  be  finite  at  r  -  0  ,  *  0 


hence 

V>  '  C1NW  *  CWU” 

The  solution  of  the  outside  annular  plate  Is: 


FN<-1  “  C1N*1F1  (uNH  ^  *  c2  N*1  F2^ +  °3  N«1F3 


CHN>-1liVUNM  ^  *  where  Vrl  “  aN-ir 
If  the  boundary  conditions  at  r  -  rNk_2  are: 


W  (r  ■  W  -  ■ "  N-2 

N**1 


—  I 

3r  1  r»r 


N^2  n’FN-2 


then  using  the  B.C.'s  at  the  mutual  Joint,  N<-1 ,  we  have: 


at  r  "  rfj“i  :  WN  N*>1  “  WN*1  N-1  “  ° 


( 


Appendix  D:  Manual  for  Buc.bas  -  buckling  analysis  of  multi-annular  plate 
(IBM  -  P.C  Version) 

Insert  In  Drive  A  the  Basic  Compiler  diskette. 

Insert  In  Drive  B  the  disk  containing  the  Buc.  exe.  and  the  input  file. 

Input  file 

Line  1:  Answer  yes  or  no  on  the  following  questions: 

1 )  Do  you  want  Intermediate  printout? 

2)  Do  you  want  normal  mode  printout? 

3)  Do  you  want  a  parametric  study?  (Two  part  only) 

4)  Is  It  buckling  in  tension  (T)  or  compression  (C)? 

Ex:  yes,  yes,  no,  T. 


Line  2: 
Case  a: 
Line  3s 
Line  4: 

yes  -  1 
Line  5: 


Input  name  of  example 

If  answer  on  question  (3)  is  ye3,  then: 

Input  lowest  and  highest  buckling  modes.  (Ex.  0,4) 

Input  to  questions:  1)  Are  there  horizontal  springs? 

2)  Are  there  vertical  springs? 
no  -  2  Ex:  2,  2,  (No  springs) 

Input  horizontal  restraints,  loading  and  springs  (If  answer  for 
question  In  line  4  Is  yes)  *  for  each  joint  on  separate  line. 


Line  5  ~  joint  0 
Line  6  -  joint  1 


Line  6  +  npr  ♦  npr: 


Input  (for  each  joint  on  separate  line) 
horizontal  restraints,  loading  and  spring  (if 
the  answer  is  yes  to  first  question  in  line  5) 

Line  6  +  npr  x  3  Input  (for  joint  on  separate  line)  vertical  b.c.'s 

and  spring  (if  the  answer  is  yes  to  second  question  in 
line  5). 

Example  (Parametric  study-  case  a) 

l  r-o  1  ane  b.  c.  -  =  =  =  =  =  =  =  =  =  =  - 

Forces:  -> 

Vertical  B.  C  -  ==========  _ 

0  1  1  0 
1  £  £  1 

1.  no,  no,  yes,  c  (no  printout,  no  printout  of  normal  mode,  parametric 
study  (two  part,  compression  force. 

2.  Example  #1  (name  of  example) 

3.  2,2  (No  springs  in  horizontal  and  vertical  directions) 

% .  0,3  (lower  mode  0,  higher  ~3) 

5.  0,1  (free  to  slide,  p  «  1  -  joint  0 

o 

6.  2,0  (Continuous  In  horizontal  direction,  p  =  0)  -  joint  1 

e1 

7.  0  B.C.  -  simply-supported  -  joint  0 

-  joint  1 

8.  10  continuous 

9.  1.5  0.3333,  1,  7  (outer  radius  -  1.5 

Poisson  ratio  -  1/3 

lowest  E2/Ej  related  to  1  and  highest  modulus 
ratios  r  7) 

\  10 


10.  1,9, 1.1  (lowest  r,/r0  -  1,  highest  Tj/Tq  -  9  and  the  factor 

multiplying  pcr  (E2/E,,  r^rQ-l)  *-1.1) 

11.  10.1,0.5,1,2,10,100  (modulus  ratios  related  to  1  to  7) 

Explanation: 

Line  7:  Input  vertical  boundary  conditions  and  spring  if  the  answer  to 
second  question  in  Line  M  is  yes. 

Line  7  -  joint  0 
Line  8  -  joint  1 
Line  9:  Input  1)  Outer  radius 

2)  Poison  ratio  (constant  for  all  parts) 

3)  Lowest  and  highest  integer  (maximum  7)  related  to  Ea/E ! 
Ex:  1.8  ,  0.3333.  2,  7 

Line  10:  Input  1)  Lowest  and  highest  integers  mintmum/maxtmum  9 
related  to  E-|/E0. 

2)  A  factor  which  multiplies  the  Pcr  found  for  a  lower 
r,/r0. 

Line  11:  Input  modulus  ratios,  E2/E,  (No.  of  values  must  be  equal  to 
(highest-lowest)  integer  appearing  in  line  9. 

Case  b:  If  answer  to  question  (3)  on  line  1  is  _no, 

Line  3:  Input  1)  No.  of  part  (No.  of  annular  +  circular  parts) 

(maximum  -  5) 

2)  No.  of  maximum  terms  in  series  (smaller  or  equal  than 


Line  4:  Input  lowest  and  highest  buckling  modes. 

Line  5:  Input  to  questions:  1)  Are  there  horizontal  springs? 


2)  Are  there  vertical  springs? 


yes,  -1 ,  no  -  2,  Ex. :  2,2 

Line  6:  Input  (for  each  joint  on  a  separate  line): 

1)  E"  modulus  of  Elasticity 

2)  P0<-  Poisson  ratio 

3)  T  -  Thickness 

Line  6  +  npr:  Input  (for  each  joint  on  a  separate  line): 
Example  (3  part  plate) 


I  x  npr 


radius  of  joint 


ir-olane  b.  c. 
Forces : 

Vertical  B.  C 


1)  no,  yes,  no,  c  (no  printout,  print  normal  modes, 

no  parametric  study,  compression) 


2) 

Example  #2 

3) 

3,  50  (no  of  partr- 

3,  no  of  terms  in  each  series  -  50 

4) 

0,5  (lowest  mode 

-  0,  highest  mode  -  5) 

5) 

2,2  (no-spring  in 

horizontal  and  vertical  direction) 

6) 

1.0.333,  1  (E(1), 

P0(1),  T(1 ))  -  plate  1 

7) 

1,0.1666,1  (E(2) , 

P0(2),  T(2 ) )  -  plate  2 

8) 

1/2,  0.0,  1  (E(3) , 

Po(3) , T ( 3 ) )  "  Plate  3 

9) 

3  (r(0) ) 

10) 

1.5  (r( 1 ) ) 

3 


:» 

\3Z 


11) 

1.0 

(r(2) ) 

12) 

0,1 

(sliding,  p  -1.0 
eo 

13) 

2,0 

(continuous  In  horizontal  dir 

HO 

2,0 

(continuous  in  horizontal  dir 

15) 

0 

(simply  supported) 

16) 

10 

(continuous) 

17)  10 


(continuous) 


Appendix  E:  Computer  Program  -Buc.bas 


10  REM  Buckling  analysis  of  circular  plate  composed  of 

annular  and  circular  part  with  different  mechanical  and 
geometrical  properties. 

20  REM  Primary  analysis  is  performed  using  stiffness 
equations  starting  at  line  1000  and  ends  at  1380 
30  INPUT  "name  of  INPUT  file  <  name  should  include  disk  drive 
letter) " ; INPU* 

40  OPEN  INPU#  FOR  INPUT  AS  #3 
50  OPEN  " 1 pt 1 : "  AS  #1 
60  REM  OPEN  0UTP$  FOR  OUTPUT  AS  #1 
70  DEFDBL  A  -H: DEFDBL  P-Z 

80  SS1$="##.  ##"  :  SS2*="###.  ###"  :  SS3$*="####.  ###"  :  IPRI*="no" 

90  DEFINT  1-0 

100  REM  Questions:  1.  Do  you  want  intermediate  printout — yes 

or  no,yepr$ 

2.  Do  you  want  normal  mode  printout — yes 
or  no,yenor$ 

3.  Do  you  want  parametric  study  (Tow 
parts) -yes  or  no,yepar* 

110  REM  "Is  it  buckling  in  "tension" (T>  or  "compression" (C) " 
120  INPUT  #3,YEPR*,YEN0R*,YEPAR*,YETY* 

130  INPUT  #3,NAMEC*:REM  INPUT  "name  of  exampl e" ; NAMEC* 

140  IF  YEPAR*=“yes"  THEN  NPR=2: NPL=NPR: NNS-80: GOTO  160 
150  INPUT  #3 , NPR , NNS : NPL-NPR : REM  input  "no.  of  elements  in 
the  power  series  solution.  No.  of  plates  including 
hoi es" ; NNS , npr 

160  INPUT  #3 , N00 , NCMAX : REM  INPUT  "The  lowest  and  highest 
buckling  mode  (in  circumf ernti al  direction" jnOO, NCMAX 
170  INPUT  #3, IYESH, IYESVsREM  INPUT  "are  there  any  spring 
support  in  horizontal  and  vertical  directions. 

Yesh=l , yesv=l , noh=2 , nov=2 " j IYESH, IYESV 
180  OPTION  BASE  0 

190  REM  Dim  declaration  for  the  primary  sub. 

200  DIM  E  (5)  ,  F'O  ( 5)  ,  T  ( 5 )  ,  BETA  ( 4>  ,  D  (5)  :  REM  starting  from  1 

210  NPLI=NPL-l:NBC=2+4*NPLIsNBC2=NBC*NBC:DIM  R (4) , IRES (4) , 

PE ( 4 ) , B ( 1 8 ) , S A (4,4) ,  S  <  4 , 4 ) ,SSL(4> ,SSR(4> , U ( 4 ) ,  AA  (324) 
,L(324> ,M(324> :  REM  -starting  from  0  , except  b  and  a 
220  IF  IYESH=1  THEN  DIM  UK (4) 

230  DIM  SOF (5) ,S0CF<5) ,SRR(5,10> ,S00(5,10) ,UU(5,10) ,RR(5,10)i 
REM  starting  from  1 

240  REM  Dim  declaration  in  lines  75-90  are  for  the  bauckling 
anal ysi s 

250  DIM  I BCW  <  4 )  ,  I BCWG* ( 4 > , WR (4,4)  , DWR (4,4) , BMR (4,4) , SHR (4,4) 
260  DIM  WL (4,4) , DWL (4,4) , BML (4,4) , SHL (4,4) 

, BMAT (18,18) , SBM (17,17) ,J00  (4) 

BMAt-starts  from  1,1  to  4#npl i +2 , 4*npl i +2 
270  DIM  PCRR (7,9) , SOFD ( 5 ) , SOCFD ( 5 ) , SRRR ( 4 ) , SRRL ( 4 ) , 

HJ ( 4 ) , JH (4) , IND (4)  ,DA(80> ,D2A(80) ,D3A<80> ,AAA(B0> , 

E22 (7) 


280  DIM  GAMMA  < 5)  ,  ALAMDA  (5)  ,  U1M  (5,  10)  ,DWW  (5,10)  ,  BMW  (5,10)  , 
SHW(5, 10)  ,  SR  (  5, 4), SMS, 4), F  (5, 4),  DF  (5, 4)  ,D2F<5,4>  , 

D3F  (5,4)  ,  D4F  <  5 , 4)  ,A(5,4,B0)  ,SU<2>  ,DSl)<2)  ,D2SU<2>  ,D3SU(2) 

, D4SU (2) , G (2) ,DG(2) ,D2G<2> ,D3G(2) ,  D4G<2) ,H(4) ,DH<4) 

, D2H (4)  , D3H (4) ,D4H (4) 

290  IF  IYESV=1  THEN  DIM  WK(4) 

300  IF  YEPAR$="yes"  THEN  GOTO  370 
310  REM  INPUT  OF  DATA 

320  REM  MECHNICAL  AND  GEOMETRICAL  PROPPERTIES 
330  FOR  J  =  1  TO  NPL: INPUT  #3 , E < J > , PO < J > , T < J > : E ( J > =E < J ) *  1 2*  ( 1- 
P0(J)*P0(J) ) 

340  D  ( J ) =E ( J  >  *T  <  J  )  "-3/  12/  <  1-PO  ( J  )  *P0  <  J  )  )  :  NEXT  J 
350  FOR  1*0  TO  NPL-ls INPUT  #3,R(I):NEXT  I 
360  FOR  1  =  1  TO  NPL-1 : BETA ( I >  =R  < 1-1 ) /R ( I ) :  NEXT  I 
370  REM  restraints  and  external  forces 

380  FOR  1=0  TO  NPL-1 J INPUT  #3,  IRES ( I ) , PE ( I ) : IF  IYESH=1  THEN 
INPUT  #3,  UK ( I ) 

390  NEXT  I 

400  REM  input  of  b.c.  for  vertical  displacement 
410  FOR  1=0  TO  NPLIs INPUT  #3 , IBCW < I > : IF  IYESV=1  THEN  INPUT 
#3 , WK ( I ) 

420  IF  I BCW  ( I )  =0  THEN  IBCWG$<I>="  A"  ELSE  IF  IBCW(I)=1  THEN 
IBCWGS  ( I )  =  "  !  "  ELSE  IF  IBCW(I>=2  THEN  IBCWG* < I > =" i i "  ELSE 
IF  IBCW  <  I )  =3  THEN  IBCWG$  ( I )  ="K-''"  ELSE  IF  IBCW(I)=4  THEN 
IBCWG4  < I ) ="  " 

430  IF  IBCW ( I ) =10  THEN  IBCWG4<I>="  "  ELSE  IF  IBCW(I>=11  THEN 
IBCWG4  ( I )  =  "  ELSE  IF  IBCW  (I )  =12  THEN  IBCWG*(I>»"  !" 

ELSE  IF  IBCW  ( I )  =13  THEN  IBCWG$  ( I )  ="K"'"  ELSE  IF  IBCW(I)  =  14 
THEN  IBCWG4 ( I ) =" ! I " 

440  NEXT  IjREM  IBCW(I>=0  to  4  means  edge  b.c.  IBCW(I)=10  to 
13  means  b.c.  at  inn  -er  joints 

450  IF  YEPAR*="no"  THEN  GOTO  530 
460  REM  "Data  for  a  parametric  study" 

470  REM  "TWO  PLATE  type  - ===== - 

480  REM  "Parameters:  l)po(poisson  rat i o) -. 3333 , . 3 , . 25 , . 1666 , 0 

2) e(l>=1.0,e(2)/e(l)=.01,.l,.5,1.2 
,10,1 00 

3) r  <0>  =  1 . 0001 ,r ( 1 > /r (0) =0. 1  to 

. 9,dec  =  .  1 

490  T ( 1 ) =1 ( : T (2) =1 !:E(1)=1 

500  INPUT  #3,R(0) ,P00, IE210, IE21F: PO < 1 ) =P00: PO <2> =P00 

510  INPUT  #3, I RIOS, IR10F,PFACT0 

520  FOR  1= IE210  TO  IE21F: INPUT  #3 , E22 ( I ) : NEXT  I 

530  FOR  NC=N00  TO  NCMAX:IF  YEPAR*="no"  THEN  GOTO  640 

540  E ( 2 ) = 1 . 1: I  YE I NTE= 1 1 GOSUB  6350: E (2) =0 

550  PRINT  " - >>>  nc="sNC: INPUT  "pO - for  nc";PCRO 

560  FOR  I E2 1  =  1 E2 1 0  TO  IE21Fi E  <2>  =E22  ( IE21 ) 

570  PRINT  "pf actO=" ; PFACTOi INPUT  "pf actO-new" ; PFACTO 
580  CCD=12* < 1-P00*P00> : D ( 1 ) =E ( 1 ) *CCD*T ( 1 ) *T ( 1 ) *T ( 1 ) / 1 2/ (1- 
P0<1)*P0<1) )  : D ( 2) =E ( 2 ) *CCD*T (2)*T (2)*T  <2) / 12/ ( 1 - 

PO (2) *PQ  <2) ) 

590  PRINT  #1 , : PRINT  #1,:PRINT  #1,:PRINT 

#1 , TAB (10) | "r 1 /rO" : TAB (25) | "e2/el=" j :PRINT  #1 , USING 
SS3*sE<2) : PR  I  NT  #1 , 


600  PRINT  "e<2>=";E<2>  ,  ,,pcr0=‘* ;  PCRO 
610  PRINT  "e  (2)  ="  ;  E  (2) ,  "pcr0=" } PCRO 
620  FOR  I R 1 0=  I R 1  OS  TO 

IR10F:R(1>=IR10/10*R(0> !BETA(1)=R(0) /R ( 1 ) 

630  PRINT  "r ( 1 ) =" j  R  < 1 ) , "pcrO=" ; PCRO 

640  REM  Main  program  which  calculates  the  Per  for  each 
nc (no.  of  circuf erenti al  waves) 

650  GOSUB  2250 

660  IF  YEPAR$="yes"  THEN  GOTO  690 

670  PRINT  "intial  value  of  F'cr=PO*DO/r  (npl  i  >  >  '2  - input  pO  - 

— for  NC=" f NCi INPUT  PO: PCR0=P0*1 ! /R (0) /R (0) : PRINT  PCRO 
680  GOSUB  1 980 : PCRR ( 0 , 0 ) =PCR 1  : GOSUB  980: GOTO  1250 
690  GOSUB  1980: PCRR ( IE21 , IR10)=PCR1: GOSUB  980: GOTO  1230 
700  FOR  1  =  1  TO  NPLI : SOFD ( I ) =PCR*SOF  < I ) /D ( I ) i IF 

ABS(SOF(I) ><=.00001  AND  ABS (SOCF ( I ) ) >. 00001  THEN 
SOF ( I >  =0: GOTO  720 

710  IF  ABS (SOCF (I) ><=.00001  AND  ABS (SOCF < I > /SOF ( I > > <=. 00001 
THEN  SOCF ( I ) =0 

720  SOCFD ( I >  =SOCF ( I>  *PCR/D ( I > : SRRR ( 1-1 > =SRR ( 1,0) *PCR: 

SRRL < I > =PCR*SRR < I  ,10): IF  SOCFD ( I > =0  THEN  PCRSI=0  ELSE 
FCONST =NC*NC- 1 

725  IF  NC=0  THEN  FC0NST2=0  ELSE  IF  NC=1  THEN  FCQNST2=-4  ELSE 
FC0NST2=-4*NC*(NC-SQR (FCONST) > 

730  IF  SOCFD  < I ) =0  THEN  PCRS I =0 : PCRS 1 2=0 : GOTO  770 
740  PCRS I =FCONST*D ( I > /SOCF ( I > : PCRSI2=FC0NST2*D ( I > /SOCF ( I ) : 
PRINT  "persi , per si  2= " a  PCRS I , PCRS 1 2: IF  PCR=PCRSI  THEN 
SOCFD ( I ) =FCONST  ELSE  IF  PCR«PCRSI2  THEN  SOCFD ( I > =FC0NST2 
750  IF  PCRSI<0  THEN  PCRSIL=PCRSI*1 . 01 : PCRSIR=. 99*PCRSI 
755  IF  PCRSI2<0  THEN  PCRSI2L»PCRSI2*1 . 01 : PCRSI2R=. 99*PCRSI2 
760  IF  PCRS I >0  THEN  PCRSIL*PCRSI*. 99: PCRSIR=1 . 01*PCRSI 
765  IF  PCRS 1 2 >0  THEN  PCRSI2L=PCRSI2*. 99: PCRSI2R=1 . 01*PCRSI2 
770  NEXT  I 

780  IF  E ( NPL ) =0  THEN  GOTO  790  ELSE  SOFD(NPL>=-  SRR(NPL,0> 

*PCR / D ( NPL  > : SOCFD  <  NPL  >  =0 : SRRR ( NPL I ) =PCR*SRR ( NPL , 0  > 

790  FOR  1=1  TO  NPL I : R I =R ( I - 1 > : IF  SOCF ( I > =0  THEN  GOSUB 
1620: GOTO  810 
800  GOSUB  3350 

810  IF  1-1=0  THEN  INDE=0  ELSE  INDE=0 
820  FOR  K= 1  TO  4: GOSUB  6130 

830  WR ( I - 1 , K ) =F (I , K > :DWR(I-1 ,K> =DW: BMR ( 1-1 ,K) =BM: SHR < I- 
1  ,K)  =SH+SRRR  (1-1)  *DW:  NEXT  K 
840  RI=R ( I > : IF  SOCF ( I > =0  THEN  GOSUB  1620: GOTO  860 
850  GOSUB  3780 

860  IF  I=NPLI  AND  E(NPL>=0  THEN  INDE=1  ELSE  1NDE=0 
870  FOR  K= 1  TO  4: GOSUB  6130 

88  •>  WL  ( I  f  K )  =F  <  I , K )  :  DWL  ( I  , K >  =DW :  BML  ( I  ,K)=BMt 
SHL  <  I  ,  K ) =SH+SRRL ( I ) *DW  : NEXT  K: NEXT  I 
890  IF  E ( NPL >  =0  OR  T(NPL>=0  THEN  GOTO  950 

900  I=NPL: RI=R (NPLI >: GOSUB  1620:F0R  KKW-1  TO  2  : IF  KKW=1  THEN 
K=1  ELSE  IF  KKW=2  THEN  K=4 
910  IF  KPROBOO  THEN  K=KKW 
920  GOSUB  6130 

930  WR ( I - 1 , KKW >  =F ( I , K  >  s DWR ( I - 1 , KKW >  =DW : BMR ( I - 1 , KKW > =BM : SHR ( I - 
1 , KKW>  =SH  +SRRR ( I - 1 >  *DW : NE  X  T  KKW 


940  REM  WR (1-1,2) =F ( I , 3 ) s  DWR (1-1,2) =DW i BMR ( I - 1 , 2  >  =BM  s SHR ( I - 
1 ,2)=SH+SRRR(I-1)*DW 
950  SOSUB  1280 
960  N=NBCC: GOSUB  5640 
970  RETURN 

980  REM  "calculation  o-f  the  Ci  -for  the  varios  •functions 
990  IF  YEPAR$="yes"  THEN  GOTO  1000  ELSE  IF  YEPAR*="no"  AND 
YEN0R*="yes"  THEN  PRINT  "nc  =  " n  NCp  TAB (30) $ 

"Per**"  ;  PCRR  (0,0)  *R  (0)  *R  (O)  /I  !  ;  "*dO/r  (O)  A2"s 
PCRR 1  =»PCRR  <  0 , 0 ) »  GOTO  1 030 

1000  BEEPsPRINT  "  nc=" « NC : PR I NT  "  e  (2)  **"  5  :  PRINT  USING 

SS2$ ; E ( 2 ) : PR I NT "  r ( 1 ) /rO=" ; s PRINT  USING 

SS2$; IR 10/ 10s  PRINT  "  Per =" 5 : PRINT  USING 

SS3$ ; PCRR ( I E2 1 , IR10) *R (0) *R (0) / 1 ! 5 SPRINT  " 

*dO/r (0) A2" s  PCRR1 =PCRR ( IE21 , IR10) 

1010  PRINT  #1 ,TAB<9> j sPRINT  #1, USING  SS2$; IR10/ 10; : PRINT 
#1 ,TAB(24> ; SPRINT  #1, USING 
SS3$ s  PCRR ( I E2 1 , IR10)*R(0)*R<0>/1 
1020  IF  YENOR*="no"  THEN  RETURN 

1030  FOR  1=1  TO  NBCC-1 : FOR  J=1  TO  NBCC-1 : SBM ( I , J > =BMAT ( I , J > s 
NEXT  J : B ( I ) =-BMAT ( I , NBCC ) s  NEXT  I 
1040  ISS=0: FOR  1=1  TO  NBCC-1: FOR  J=1  TO  NBCC-1 : ISS=ISS+1 : 

AA ( ISS) =SBM ( J , I ) : NEXT  JsNEXT  I s N=NBCC- 1 s GOSUB  2910s 
B (NBCC) =1 ! s GOTO  1050 

1050  FOR  1=1  TO  NBCCs PRINT  "b ( " ; I ; " ) =" ; B ( I ) « NEXT  I 
1060  FOR  1=1  TO  NPLI : FOR  JW=0  TO  10s RI =RR ( I , JW) : IF  1=1  AND 
JW=0  THEN  INDE=1  ELSE  INDE«0 
1070  IF  I=NPLI  AND  JW=10  AND  E(NPL)=0  THEN  INDE=1  ELSE 
INDE-0 

1080  IF  SOCF ( I ) =0  THEN  GOSUB  1620  ELSE  GOSUB  3350 
1090  INJ=4* ( I  — 1 ) : SUMW=0:  SUMDW=0: SUMBM=0: SUMSH=Os  FOR  K=1  TO 
4  s  I N J  K= I N J +K  3  GOSUB 

6130s  SUMW=SUMW-*-B ( INJK) *F ( I ,K) s SUMDW=SUMDW+B < INJK) *DF ( I , K) s  SUM 
BM=SU  MBM+B< INJK) *BM: SUMSH-SUMSH+B ( INJK) *SH: NEXT  K 
1 1 00  WW ( I , J W ) =SUMW  s  DWW  < I , J W ) -SUMDW  t BMW ( I , J W ) =SUMBM : 

SHW ( I , JW) =SUMSH  +SRR (  I , JW) *PCRRl*SUMDWs NEXT  JWsNEXT  I 
1110  IF  E  < NPL ) =0  THEN  GOTO  1160 

1120  I=NPL: FOR  JW=0  TO  5  s RI=RR < I , JW) s IF  JW=5  THEN  RI*. 00001 
1130  IF  SOCF ( I ) =0  THEN  GOSUB  1620  ELSE  GOSUB  3350 
1140  INJ=4*NPLI : SUMW=0: SUMDW=0: SUMBM=0: SUMSH=0: FOR  K=1  TO  4 
STEP  3s INJ=INJ+1 s GOSUB 

6130s  SUMW=SUMW+B  <  I N J )  *F  ( I  ,  K )  s  SUMDW=SUMDW-t-B  ( I N J  )  *DF  ( I  ,  K )  s  SUMBM 
=SUMB  M+B ( INJ ) *BMs  SUMSH=SUMSH+B ( INJ ) *SHs  NEXT  K 
1 150  WW( I , JW) =SUMWs DWW ( I , JW) =SUMDWs BMW ( I , JW)=SUMBMs 
SHW(I ,JW)=SUMSH  +SRR<  I , JW) *PCRRl*SUMDWs NEXT  JW 
1160  WMAX=0: FOR  1=1  TO  NPLIsFOR  J=1  TO  10: IF  ABS ( WW ( I , J ) ) >= 
ABS(WMAX)  THEN  WMAX=WW(I,J) 

1170  NEXT  JsNEXT  I 

11 SO  IF  E (NPL) =0  THEN  GOTO  1210 

1190  I=NPLs  FOR  J  =  1  TO  5s  IF  ABS ( WW ( I , J ) ) >=ABS ( WMAX >  THEN 
WMAX=WW  < I , J) 

1200  NEXT  J 
1210  GOSUB  7030 
1220  RETURN 


1230  PCR0-PFACT0*PCRR(IE21 , I  RIO) sNEXT  IR10 
1240  PCRO*PCRR ( IE21 , IR10S) : NEXT  IE21sG0SUB  6350 
1250  NEXT  NC 
1260  CLOSE 
1270  END 

1280  REM  Sub.  which  constructs  the  Buckling  Matrix 

1290  FOR  1*0  TO  NPLIsIF  1*0  THEN  GOTO  1450 

1300  IR0W=2+ ( 1-1 ) *4: ICOLM* ( I- 1 ) *4: Il*IROW+ls I2»IRQW+2i 

1 3* I RQW+3 c 1 4= 1 ROW  +4: IF  I*NPLI  THEN  KMAX*6  ELSE  KMAX-8 
1310  IF  I=NPLI  AND  E(NPL)=0  THEN  GOTO  1510 
1320  IF  I BCW ( 1 ) * 1 0  THEN  GOTO  1330  ELSE  IF  IBCW(I)*11  THEN 
GOTO  1360  ELSE  IF  IBCW(I>=12  THEN  GOTO  1390  ELSE  IF 
I BCW  <  I )  *13  THEN  GOTO  1420 

1330  FOR  K=1  TO  KMAXs ICK=ICOLM+Ks IF  K>4  THEN  GOTO  1340  ELSE 
BMAT (II, ICK) =WL ( I ,K) s BMAT ( 12, ICK) *DWL  < I ,K>  » 

BMAT (13, ICK) *BML ( I  ,K)sB  MAT ( 14, ICK) =SHL ( I ,K> : GOTO  1350 
1340  K4»K-4j  BMAT  <11, ICK) *-WR ( I ,K4) s  BMAT ( 12, ICK) *-DWR ( I ,K4) s 
BMAT ( 1 3 , I CK )  *-BMR ( I ,K4) : BMAT  < 1 4 , 1 CK ) =-SHR ( I , K4 ) 

1350  NEXT  Ks GOTO  1570 

1360  FOR  K-l  TO  KMAXj ICK=ICOLM+Ks IF  K>4  THEN  GOTO  1370  ELSE 
BMAT (II, ICK) =WL  < I ,K> : BMAT ( 12, ICK)*0 ! s 
BMAT  <13, ICK) =DWL ( I ,K)  sBMAT(I4  , ICK) *BML ( I ,K> I 
GOTO  1380 

1 370  K4*K-4 : BMAT (11,1 CK ) *0 ! i BMAT ( 1 2 , I CK >  =WR ( I ,K4) : 

BMAT ( I 3 , I CK ) =-DWR ( I , K  4) s BMAT ( 14 , ICK) *-BMR ( I ,K4) 

1380  NEXT  KjGOTO  1570 

1390  FOR  K=1  TO  KMAXs ICK*ICOLM+Kl IF  K>4  THEN  GOTO  1400  ELSE 
BMAT (II, I CK ) *WL  < I , K ) s BMAT ( 12 , ICK) *0 ! s 
BMAT ( 1 3 , I CK ) *DWL  < I , K ) s  BM  AT (14  , ICK) *0 ! s GOTO  1410 
1400  K4=K-4s  BMAT ( 1 1 , ICK) *0 ! : BMAT ( 12, ICK) *WR ( I ,K4) s 
BMAT ( 13 , I CK >  *0 ! s BMAT (  14, ICK) »DWR < I ,K4> 

1410  NEXT  K: GOTO  1570 

1420  FOR  K* 1  TO  KMAX: I CK* I COLM+K : IF  K>4  THEN  GOTO  1430  ELSE 
BMAT (11,1 CK ) *WL ( I , K ) s BMAT ( 12, ICK) *DWL ( I ,K) s 
BMAT ( 13, ICK) *BML ( I  ,K):B  MAT < 14, ICK) *SHL < I ,K> - 
WK ( I ) *WL  < I , K ) : GOTO  1440 

1 430  K4*K-4 : BMAT ( 1 1 , ICK ) *-WR ( I ,K4) i BMAT ( 1 2 , I CK ) =-DWR ( I ,K4) s 
BMAT (13, ICK)  *-BMR ( I , K4 ) s BMAT ( 14, ICK) »-SHR ( I ,K4) 

1440  NEXT  K: GOTO  1570 

1450  IF  I BCW (0) =0  THEN  GOTO  1460  ELSE  IF  IBCW(0)*1  THEN  GOTO 
1470  ELSE  IF  IBCW(0)*2  THEN  GOTO  1480  ELSE  IF  IBCW(0)-3  THEN 
GOTO  1490  ELSE  IF  IBCW(0)=4  THEN  GOTO  1500 
1460  FOR  K= 1  TO  4s BMAT ( 1 , K) =WR <0, K) s BMAT (2 , K> *BMR (0,K> l NEXT 
K s GOTO  1570 

1470  FOR  K* 1  TO  4s BMAT ( 1 , K ) =WR (0 , K) i BMAT (2 , K) =DWR (0, K) s NEXT 
Ks GOTO  1570 

1480  FOR  K=1  TO  4s BMAT ( 1 , K > =DWR (0 , K) » BMAT <2 , K) *SHR <0 ,K> s NEXT 
K s GOTO  1570 

1490  FOR  K=1  TO  4s BMAT ( 1 , K) =BMR <0, K) s BMAT < 2, K> =SHR <0 ,K) - 
WK ( I ) *WR (0 , K) s  NEXT  KsGOTO  1570 
1500  FOR  K-l  TO  4s BMAT ( 1 ,K) =BMR <0,K> s BMAT (2,K> *SHR <0,K> s NEXT 
KsGOTO  1570 


1510  IF  IBCW  ( I )  =1 1  THEN  GOTO  1520  ELSE  IF  IBCW(I>=12  THEN 
GOTO  1530  ELSE  IF  IBCW(I>“13  THEN  GOTO  1550  ELSE  IF 
IBCW ( I ) *14  THEN  GOTO  1540  ELSE  IF  IBCW(I)=10  THEN 
GOTO  1560 

1520  FOR  K-l  TO  4«  BMAT  (II,  ICOLM+K)  =WL  ( I , K) l 

BMAT ( 12, ICOLM+K) “BML ( I ,K) i NEXT  KiGOTO  1570 
1530  FOR  K“1  TO  4» BMAT < 1 1 , ICOLM+K) »WL ( I ,K> s 

BMAT ( 12, ICOLM+K) “DWL  < I ,K) s  NEXT  KsGOTO  1570 
1540  FOR  K-l  TO  4j BMAT ( 1 1 , ICOLM+K) -DWL < I ,K> ! 

BMAT <12, ICOLM+K )=SHL< I, K) sNEXT  KiGOTO  1570 
1550  FOR  K»1  TO  4s BMAT < 1 1 , ICOLM+K) -BML ( I ,K> s 

BMAT  <12,1 COLM+K ) =SHL  < I , K ) -WK ( I ) *WL ( I , K  ):NEXT  K: 

GOTO  1570 

1560  FOR  K=1  TO  4s BMAT ( 1 1 , ICOLM+K) »BML ( I ,K> s 

BMAT (12, ICOLM+K) =SHL< I, K) sNEXT  KsGOTO  1570 
1570  NEXT  I 

1580  IF  E (NRL) *0  THEN  NBCC=4*NPLI  ELSE  NBCC=2+4*NPLI 
1590  ISS*Os  FOR  1*1  TO  NBCCsFOR  J  =  1  TO  NBCCs ISS=ISS+1 : 

AA ( ISS) “BMAT ( J , I ) s  NEXT  JsNEXT  I 
1600  IF  YEPR*="no"  THEN  GOTO  1610  ELSE  FOR  1=1  TO  NBCCsFOR 
J=1  TO  NBCCs PRINT  "bmat ( " 5  I ; " , " s J " > 5 BMAT ( I , J > 

, ; SNEXT  J : PRINT : NEXT  I 
1610  RETURN 

1620  REM  SUB.  "solves  the  D.E.  for  a  uniform  state  of  in¬ 
plane  stresses"-  srr=-sof ( i > *pcrr ( i > 

1630  KPR0B=0i SR ( 1 , 1 ) -2+NCs  SR < 1 ,2) -2-NCs  SR ( 1 ,3) =-NCs 

SR (1,4) *NC  sSI(I,l)*0  ! SI ( 1 , 2) *0s  SI ( I , 3) *0s  SI (1,4) =0 
1640  PRINT} TAB (10) }"  i  ";TAB(30>;"  sr  "sTAB(45)j"  si 
"SPRINT 

1650  FOR  K=1  TO  4s  PRINT; TAB ( 12) s 1 1 TAB (30) s SR < I , K) 

; TAB (45) ; SI (I ,K) sNEXT  Ks  PRINT 
1 660  K* 1 s  JUP=NNS*2 : NNS2“2*NNS 
1670  A(I,K,0)=ls  FOR  J*2  TO  JUP  STEP  2sJ2“J/2 
1680  DEN0M2=J+SR< I ,K) s DEN0M»DEN0M2*DEN0M2-NC*NC 
1690  A  <  I , K  ,  J2)  —SOFD  ( I )  *1  /DEN0M*A  ( I  ,K  ,  J2-1 ) 

1700  REM  PRINT  " i =" ; I , " k= " s K , " j=" s J , "a ( i , k , j ) = " s A ( I , K , J2) 
1710  NEXT  J s IF  K=2  THEN  RETURN 
1720  KPP=Ks  JO=Os  GOSUB 

4730s  F ( I , K ) “SUM s DF ( I ,K> =DSUMs  D2F ( I , K > =D2SUMs  D3F ( I . K) =D3SUMs  D4 
F  < I , K  ) =D4SUMs  GOSUB  6290s  IF  YEPR*“"yes"  THEN  PRINT 
"F(I,K),K  " ; F ( I , K ) , K 

1730  K“2s  IF  NOO  THEN  GOTO  1770  ELSE  EPSS-.  0000001 : 

SR ( 1 ,2) “SR (1,1 )+. 0000001 s  JUP=NNS2s  A < I ,K,0) =1 ! s 
GOSUB  1670 

1740  KPP*Ks  JO=Os JMAX=NNSs  GOSUB  4730 

1 750  FE1 =SUMs  DFE1 -DSUMs D2FEl=D2SUMs D3FEl“D3SUMs  D4FE1 “D4SUM 
1760  F  < I ,2)  =  (FEl-F (1,1)) /EPSSs  DF ( I ,2) “ (DFE1-DF (1,1) ) /EPSSs 
D2F ( I ,2)  =  (D2F  E1-D2F (1,1)) /EPSSs  D3F ( I ,2)  =  (D3FE1- 
D3F ( I , 1) ) /EPSSs  D4F ( I , 2) “ (D4  FE1-D4F ( I , 1 ) > /EPSS 
s GOSUB  6290s  IF  YEPR*="no"  THEN  GOTO  1910  ELSE 
PR I NT "f (i , k ) ,k"|F(I,K) , KsGOTO  1910 
1770  K“2s  IF  NC“  1  THEN  C2-1  !  s  Cl—  C2*S0FD  ( I  > /A  ( I  ,  1 ,0) /4/  ( 1+NC) 
|A(I,K,1)=1 ! i A  < I , K , 0) “Os  GOTO  1820 
1780  K“2s IF  NC=2  THEN  A ( I , K , 0) =1 ! s A < I , K , 1 > “Os C2=-S0FD ( I > /4s 
A(I,K,2)*1 ! sCl“-C2*S0FD(I>/4/  ( 1+NC> s GOTO  1820 


1790  K>2i  JUP*2*NC-42  GOSUB  1670 

1  BOO  J«=2*NC-2 1  J2=J/2i  C2*A  ( I ,K, J2-1 )  *SQFD  <  I  >  *  <  ( NC-2 )  *  ( NC-2 )  - 
NC*NC)/NC/(N  C-l) /82 A(I ,K,J2>*0 
1810  J=2*NC: J2=J/2sCl=-C2*S0FD< I ) /4/ ( l+NC) /A (I , 1 ,0) s 
A  ( I  ,  K  ,  J2>  *  1  ! 

1820  FOR  J*2*NC+2  TO  NNS2  STEP  2: J 1=J-NC*2: J J2=J 1 /2: J2=J/2 
1830  DAA1*A(I,1,JJ2)*(J1+SR(I,1)-1)*(4*(JH-SR(I,1> 

) *  <  J 1 +SR (1,1) -2) -4*NC  *NC-SOCFD  < I >  *2) 

1840  DAA2=2*S0FD  < I ) * ( J 1+SR < 1 , 1 ) -2 >  *A ( I ,  1 , J J2-1 ) 2 
DA 1 =D  AA 1 +DAA2 

1850  ANOM* (J-NC) * (J-NC) -NC*NC 2  DENOM* ( J  +2-NC ) * ( J +2-NC > - 
NC*NC s  A ( I , K , J2 )  —  ( C 1 *DA 1 / ANOM+A ( I ,K, J2-1 ) *SOFD ( I ) > 

/DENOM2  NEXT  J 

1860  SMUL*1 2  FOR  J*1  TO  NC2  SMUL*SMUL*RI 1  NEXT  J 
1870  F  ( 1 , 3 ) *C 1 *F (1,1) +C2*SMUL  2 

DF ( I ,3) »C1*DF (1,1) +C2*NC*SMUL/RI 2 
D2F (1,3) *C1*D2F (1,1) +SMUL*NC* (NC-1 ) *C2/RI /RI 
1880  D3F (1,3) *C1*D3F (1,1) +C2*NC* (NC-1 >  * (NC-2) *SMUL/R I /RI /RI : 
D4F (1,3) *D4  F ( I , 1 >  *C1+C2*NC* (NC-1 >  * (NC-2) * (NC-3) * 

SMUL/RI /RI/RI/RI 

1890  K'*22  JMAX*NNS:KPP=K2 GOSUB  4730 
1900  K=2  2  KKK=3 1  I  DS=  1 2  GOSUB  5050 

1910  IF  R 1*0  THEN  GOTO  1930  ELSE  IF  NC*0  THEN  GOTO  1920  ELSE 
K*3  2  POW=-NC : GOSUB  19502  GOTO  1940 
1 920  I DS= 1 2  SUM«0  2  D2SUM-0  2  DSUM*0 : D3SUM=0  2  D4SUM-0  2 

K-3  2  KKK=4  2  F  ( I ,  KKK ) = 1 2  DF  ( I ,  KKK )  *0 ;  D2F  ( I ,  KKK )  *0  2 
D3F  ( I ,  KKK)  =02  D4F  ( I  ,  KK.K)  *02  GOSUB  5050  2  GOTO  1940 
1930  K>3  2  F ( I , K ) "0 1 DF  < I , K ) *0  2  D2F ( I , K ) *0  2  D3F (I , K ) *0  2 

D4F ( I , K ) *0 2  GOSUB  62902  IF  YEPR$»"no"  THEN  GOTO  1940  ELSE 
PRINT "f (i ,  k) ,k"»F(I,K) ,K 
1940  K*4 2  POW*NC 2  GOSUB  1950: RETURN 

1 950  F ( I ,K)*RI ''■POM 2  DF  ( I ,  K )  «POW*R I  '  ( POW- 1 )  1  D2F  ( I ,  K)  *POW*  ( POW- 
1 )  *RI'N  (POW-2)  2  D3F  ( I  ,  K)  “POW*  (POW-1 )  *  (POW-2)  *RI^  (POW-3) 

1 960  D4F ( I , K ) *POW* ( POW- 1 ) * ( POW-2 ) * ( POW-3 ) *R I  '  ( POW-4  > 

1970  GOSUB  62902  IF  YEPR*="no"  THEN  RETURN  ELSE  PRINT" f  ( i , k ),  k 
" ; F ( I ,K) ,K2 RETURN 

1980  REM  Sub.  which  calculates  the  Per  using  Newton-Raphson 
iterative  algoritmus 

1990  EPS*. 00001 2 IEND=302  IER=02  DPCR=. 00001 2 ISING=0 
2000  PCR*PCRO: GOSUB  7002 DET0=DET2 BEEP2 PR INT "det , per " , DET , PCR 
2010  IF  ABS (DET) <*.000001  THEN  PCR1=PCR2 RETURN 
2020  FOR  1 1 IT*1  TO  I END 

2030  IF  PCRO>*PCRSIL  AND  PCR0<=PCRSIR  THEN  I S I NG* 1  + 1 S I NG  2  GOTO 
2180 

2035  IF  PCR0>*PCRSI2L  AND  PCR0<=PCRSI2R  THEN 
ISING*1+ISING2G0T0  2235 

2040  PCR*PCR0+DPCR2G0SUB  7002  DDET*  (DET-DETO) /DPCR.- IF  DDET*0 
THEN  GOTO  2170  ELSE  PCR1=PCR0-DET0/DDET 
2050  IF  PCR 1 < 0  THEN  PCRl*PCR0/2 
2060  PCR*PCR 1 2  GOSUB 

700 2  DET 1 =DET 2  BEEP 1  PR  I NT " pc r 1 , det 1 " , PCR 1 , DET 1 

2070  IF  ( ABS ( (PCR1 ) ) >1  AND  ABS ( (PCR1-PCR0) /PCR1 ) <*100*EPS> 

AND  ABS  ( DET  IX  *EPS  THEN  RETURN 


2080 


2090 


2100 


2110 


2120 

2130 

2140 

2150 

2160 

2170 


2190 

2200 

2210 

2220 

2230 


2235 


2240 

2250 

2260 

2270 

2280 

2290 

2300 


2320 

2330 

2340 

2350 

2360 

2370 


IF  (ABS (PCR1 >  <  =  1  AND  ABS (PCR1-PCR0) <»100*EPS>  AND 
ABS ( DET 1 ) <  =EPS  THEN  RETURN 

IF  ( ABS  ( PCR 1 )  <  =ABS  ( PCRS I )  AND  ABS  (F'CRO)  >=ABS  (PCRSI  >  )  OR 
(ABS  <PCR1 ) >®ABS (PCRSI )  AND  ABS (PCRO) C-ABS (PCRSI > )  THEN 
ISING*ISING+1 : GOTO  2180 

IF  DET1*DET0<0  THEN  PCR* <DET1*PCR0-DET0*PCR1 > / (DET1- 
DETO) : GOSUB  700: PCR1=PCR: DET1=DET: IF  ABS (DET) <-. 000001 
THEN  RETURN 

IF  ( ABS ( PCR 1 ) >= . 95*ABS ( PCR20 )  AND 

ABS  < PCR 1 ) <  =  1 . 05*ABS ( PCR20 ) )  AND  ABS (DET1 ) >-. 01  THEN 
PCR 1 « ( PCR 1 +PCRO ) / 2 1 PCR=PCR 1 : GOSUB  700 s  DET 1 =DET s  BEEP  s 
PRINT  "Bi-section  meth . -per  1 , det 1  ":PCR1,DET1 
PCR20*PCR0s PCR0*PCR1 : DET0-DET1 

BEEP: BEEP: PRINT"pcrO=" ; PCRO, "detO*" ; DETO , "ii it*"; IIIT: IF 
ABS  ( PCRO  )  <  * .  05  THEN  GOTO  2240 

IF  (ABS ( PCR 1 -PCRO) <=.01  AND  ABS (DET1 ) <-. 01 )  AND 
1 1  IT* I END  THEN  RETURN 
NEXT  IIIT 

PRINT  "there  is  no  converqance  after  "5IEND5" 
loops"! GOTO  2240 

PRINT  "the  derivative  is  equal  to  <<<<  zero  >>>  .  Try 

another  intial  Per": GOTO  2240 

PCRO=PCRS I : PCR=PCRO  s  GOSUB  700 : DETO* DET : BEEP  s 

RRINT"detsi ,pcrsi ", DET, PCR 

IF  ABS (DET) 0.000001  THEN  PCR 1*PCR: RETURN 

IF  ISING-1  THEN  GOTO  2230 

IF  ABS  (PCRO)  >*1  THEN  PCR0*3*PCRSI  ELSE  IF  PCR8K0  THEN 
PCRO—  ,53+PCRO  ELSE  PCRO*. 5S+PCR0 
GOTO  2000 

BEEP: BEEPs BEEP: PRINT  "Per  is  in  PCRSI*" ; PCRSI ; " 

neighborhood - <<<<<  please  enter  your  guess 

- at  least  far  from  persi - ": INPUT  PCRO: GOTO  1980 

BEEP: BEEP: BEEP: PRINT  "Per  is  in  PCRS I 2= " ; PCRS 12;" 
neighborhood - <<<<<  please  enter  your  guess 


1 : INPUT  PCRO: GOTO  1980 
guess  for  the  case 


using  stiffness  equations 


BEEP: BEEP: BEEP: PRINT  "Per  is  in  PCRSI*" ; PCRSI ; " 

neighborhood - <<<<<  please  enter  your  guess 

- at  least  far  from  persi - ": INPUT  PCRO: GOTO  1980 

BEEP: BEEP: BEEP: PRINT  "Per  is  in  PCRS 1 2=  " : PCRS 12;" 

neighborhood - <<<<<  please  enter  your  guess 

- at  least  far  from  persi - ": INPUT  PCROsGOTO  1980 

INPUT  "<<<<<  input  your  intial  guess  for  the  case 
>>>>"? PCRO: GOTO  1980 

REM  primary  analysis  using  stiffness  equations 
REM  stiffness  matrix 
FOR  1*0  TO  NPL-1 
IF  I -NPL-1  GOTO  2320 

SSR ( I ) =E ( 1  +  1 ) *BETA  < 1  +  1 >  * ( 1+PO ( I  +  l )  +  ( 1-PQ ( 1  +  1 ) > 

/BETA  ( I  + 1 )  ''2 >  /R  <  I  + 1  >  /  <  BETA  <  I  + 1  >  '  2- 1 )  /  ( 1  -PO  ( I  + 1 )  ~2 ) 

S(I  +  1,  I)-E(I  +  1)*BETA(I  +  1)*2/R(I  +  1>  /  (BETA  (  H-l  >  ~2-l  >  /  <  1  - 
PO ( 1+1 ) A2) :  SA<Z+1 , I)=S(I+1 , I> 

IF  1*0  THEN  S (0 , 0) *-SSR ( I ) : GOTO  2370 
SSL  ( I )  =E  ( I  >  *  ( 1 +P0  ( I )  +BETA  ( I )  y’2*  ( 1  - 
PO  ( I )  )  )  /R  ( I )  /  (BETA  ( I )  ''2-1 )  /  (  i-PO  ( I )  "2) 

S  <  I  — 1 ,  I )  =E  ( I )  *2/R  ( I )  /  ( BETA  ( I  )  •  "2- 1 )  /  ( 1  -PO  (I)  -s2):SA(I- 
1 , I >  =S  < I -1 , I > 

IF  I=NPL-1  GOTO  2350  ELSE  GOTO  2360 
SSR (NPL-1 ) *E (NPL  > /R (NPL-1 ) / ( l-PO(NPL) ) 

S(I  ,  I)=-SSR(I)-SSL(I) 

SA(I , I)*S( I , I) i NEXT  I 


2380  FOR  1*0  TO  NPL-1 : B < 1  +  1 > -PE ( I ) : NEXT  I 
2390  FOR  J  =0  TO  NPL-1 

2400  IF  I RES ( J ) =0  OR  IRES(J>=2  THEN  GOTO  2430 
2410  FOR  1*0  TO  NPL-1 1  SA ( I , J' *0: SA (J , I ) =0: NEXT  I 
2420  SA ( J , J ) =S (J,J) :  B  ( J  + 1 ) *0 
2430  NEXT  J 

2440  ISS*Os  FOR  1*0  TO  NPL-lsFOR  J»0  TO  NPL-1 : ISS* ISS+1 : 

AA  ( ISS) =SA  < J , I ) : NEXT  JsNEXT  I 
2450  REM  solution  of  the  deformation 
2460  N*NF'L:  60SUB  2910 

2470  FOR  1*0  TO  NPL-1 : U < I ) *B ( 1+1 ) s NEXT  I 

2480  IF  YEPR*="yes"  THEN  GOSUE'  3160 

2490  REM  calculation  of  stresses  and  deformation 

2500  FOR  1*1  TO  NPL-1 

2510  SOF  < I ) *E ( I >  * ( 1 +P0 ( I ) >  * (U ( I > -BETA ( I ) *U ( 1-1 ) >/R<I) 

/  ( BETA  <  I )  "-2- 1  >  /  <  1  -PO  ( I  >  •'■•2  > 

2520  SOCF ( I ) *BET A < I ) *E  < I ) *  < 1 -PO ( I ) >  * ( -BETA (I)*U(I)+U(I- 
1 >  >  *R  ( I  >  /  (BETA  ( I  )  --2-1  >  /  ( 1-F'O  ( I )  '-2) 

2530  DEC* ( R ( I - 1 ) -R ( I ) ) / 1 0 

2540  FOR  J*0  TO  10 

2550  RR ( I , J  >  *R  < I - 1 > -DEC* J 

2560  SRR ( I , J ) *-SOF ( I ) +SOCF ( I ) /RR ( I , J ) /RR ( I  ,  J  > 

2570  SOO ( I , J ) =-SOF ( I ) -SOCF ( I ) /RR ( I , J ) /RR ( I , J ) 

2580  UU  ( I  ,  J  )  -  <U  ( 1-1 )  *BETA  <  I  >  *  ( 1-R  ( I )  -'2/RR  ( I  ,  J  )  "'2)  +U  ( I  >  *  (R  ( I- 
1 )  "'2/RR  ( I ,  J  )  '"2-1 )  )  *RR  ( I ,  J)  /R  <  I )  /  (BETA  ( I  >  ■"'2-1) 

2590  IF  I>0  AND  J*0  AND  ABS (SRR ( I ,0) -SRR ( 1-1 , 10) ><*. 0000001# 
THEN  SRR (1,0) *SRR (1-1,10) 

2600  NEXT  JsNEXT  I 

2610  REM  stresses  and  deformation  in  circular  part 
2620  IF  E ( NPL ) =0  THEN  GOTO  2690 
2630  FOR  J=0  TO  5 

2640  SRR (NPL , J ) *S (NPL) *U (NPL-1 ) /R (NPL- 1 ) / ( 1-PO (NPL) ) : IF  J=0 
AND  ABS < SRR < NPL, J) -SRR (NPL I, 10) ><*.0000001#  THEN 
SRR ( NPL , 0 ) =SRR ( NPL I , 10) 

2650  SOO ( NPL , J ) =SRR ( NPL , J ) 

2660  RR ( NPL , J ) =R ( N- 1 ) -R ( N- 1 ) * J / 5 
2670  UU ( NPL , J ) =U  <  N- 1 ) / R ( I - 1 ) *RR ( NPL , J ) 

2680  NEXT  J 


2690  REM  printout  of  stresses 

2700  IF  YEPRf ="no"  THEN  RETURN 

2710  IF  E ( NPL ) =0  THEN  NPLMAX=NPL-1  ELSE  NPLMAX=NPL 

2720  FOR  1*1  TO  NPLMAX 

2730  PRINTsPRINT  "Plate  no.  :  " $  I s PRINT: PRINT: PRINT  : IF 

I=NPL  GOTO  2760 

2740  SSP$="###. ###" : PRINT  TAB (20) ; "no ( " ; I s " ) =  " s : PRINT  USING 
SSP$; SOF (I) SPRINT  TAB (20) ; "noc ( " s 1 5 " ) =" s s PRINT  USING 
SSP*:SOCF(I) 

2750  IF  SOF ( I ) =0  THEN  GOTO  2760  ELSE  PRINT  TAB (20) 5 

"betac  (";  I;  ")*"  s  :  PRINT  USING  SSF't ;  SOCF  ( I )  /  SOF  ( I ) 

2760  IF  I  =NF'L  THEN  SJ*5  ELSE  SJ=10 

2770  PRINT  sF'RINT  SPRINT  TAB  (10)  5"  rr  "  ;  TAB  (25)  5"  nrr 
” ;  TAB  ( 45)  ;  "  noo  "»TAB(63)s"  uu  "  s  PRINT.-  PR  INT 

2780  FOR  J*0  TO  SJ 


2790  SSP#="###. ###" : PRINT  TAB (9 ) s : PR INT  USING  SSP* | RR ( I , J > 5  s 
PRINT  TAB (25) ; : PRINT  USING  SSP*; SRR ( I , J > » t PRINT 
TAB (45) s SPRINT  USING  SSP$; S00 < I , J ) ; : PRINT  TAB (62) 5 : PRINT 
USING  SSP$;UU(I ,J) 

2800  NEXT  J : NEXT  I 
2810  FOR  1=1  TO  NPLMAX 

2820  PRINT  #1, SPRINT  #1, "PI  ate  no.  :  "jIsPRINT 

#1, SPRINT  #1, SPRINT  #l,sIF  I=NPL  GOTO  2850 
2830  SSP$="###. ###" s  PRINT  #1,  TAB (20) ; "no ( " ; I ; " ) =" 5 s PRINT  #1, 
USING  SSP$; SOF ( I ) s  PRINT  #1,  TAB (20) ; "noc ( " 5  1 5 " ) =" 

; SPRINT  #1,  USING  SSP$«SOCF(I) 

2840  IF  SOF ( I ) =0  THEN  GOTO  2850  ELSE  PRINT  #1, TAB (20); 

"betac ("; I; ")="; SPRINT  #1,  USING  SSP* ; SOCF ( I ) /SOF ( I > 

2850  IF  I=NPL  THEN  SJ=5  ELSE  SJ=10 

2860  PRINT  #1, sPRINT  #1, sPRINT  #1,  TAB (10);"  rr  " ; TAB (25) ; " 
nrr  ";TAB(45);M  noo  " ; TAB (63) s "  uu  "sPRINT  #1, 
SPRINT  #1 , 

2870  FOR  J-0  TO  SJ 

2880  SSP$="###. ###" 5  PRINT  #1,  TAB (9) ; s PRINT  #1,  USING 
SSP*;RR ( I ,J) ; sPRINT  #1,  TAB (25) ; s PRINT  #1,  USING 
SSP*sSRR(I,J) 5 sPRINT  #1,  TAB (45) j s PRINT  #1,  USING 
SSP$; SOQ ( I , J ) s  s  PRINT  #1,  TAB (62) ; s PRINT  #1,  USING 
SSP#;UU(I ,J) 

2890  NEXT  JsNEXT  I 
2900  RETURN 

2910  REM  sub.  -for  solving  simulteneous  equati  ons  (SIMO) 

2920  T0L=0  s  KS=0  s  J  J =-N 

2930  FOR  J=1  TO  Ns JY=J+1 s J J=J J+N+l : BIGA=Os IT=J J-J 
2940  FOR  I=J  TO  Ns  I J  =  IT+1 5 BA-ABS (BIGA) -ABS (AA  ( I J ) ) 

2950  IF  8A< 0  THEN  GOTO  2960  ELSE  GOTO  2970 
2960  BIGA=AA ( I J ) s I MAX = I 
2970  NEXT  I 

2980  BA=ABS (BIGA) -TOLs IF  BA<=0  THEN  GOTO  2990  ELSE  GOTO  3000 

2990  KS= Is  RETURN 

3000  Ti=J+N*<J-2) . IT=IMAX-J 

3010  FOR  K=J  TO  Ns  1 1  =  1 1+Ns 12=1 1  +  ITs SAFE=AA < I  1 ) s 
AA (II) =AA (12) s  AA ( 12) =SAFE 
3020  A A (II) =AA (II) /BIGAs  NEXT  K 

3030  SAFE=B ( IMAX ) s B ( IMAX ) =B ( J ) ; B ( J ) =SAFE/ BIGAs  JN=J-N 
3040  IF  JNOO  THEN  GOTO  3050  ELSE  GOTO  3110 
3050  IQS=N*(J-1) 

3060  FOR  I X=JY  TO  Ns  I  X J= IOS+ 1  X s I T=J- 1 X 
3070  FOR  J X=JY  TO  Ns  I  X J X=N* < J X-l ) + 1 X s J J X= I  X J X  + 1 T s 
AA ( I X J X ) =AA (IXJX)-AA ( I  X J ) *AA  < J J X ) 

3080  NEXT  JX 

3090  B ( I X  >  =B ( I X  > -B  < J ) *AA  < I  X J ) s  NEXT  IX 

3100  NEXT  J 

3110  NY=N~ 1 s IT=N*N 

3120  FOR  J=1  TO  NYs IA=IT-Js IB=N-Js IC=N 

3130  FOR  K=1  TO  Js B < IB) =B ( IB) -AA ( IA) *B ( IC> 5 IA=IA-Ns IC=IC- 
1 s  NEXT  Y 
3140  NEXT  J 

3150  RETURNS  'END  OF  SUBROUTINE 
3160  REM  printout  o-f  result 


3170  I PR  I " yes " : GOSUB  6390 

3180  PRINT: PRINT: PRINT"  <<<<<<<<<  data  and  results 

>>>>>>>> >>>>"! PRINT 

3190  PRINT  # 1 , : PR  I NT  # 1 , : PR I NT  # 1 , "  <  <  <  <  <  <  <  <  <  data  and 
results  >  >  >  >  >  >  >  >  > > >  > " : PR I NT  # 1 , 

3200  PRINT  TAB (5) ; "plate  no.  ” ; TAB ( 18) ; "  e  "  ;  TAB (30) ; "  po 
" ; TAB  ( 45) ; "  th  ";TAB<62>;"  beta  "sPRINT 
3210  PRINT  #1,  TAB (5) ; "plate  no.  " j TAB ( 18) ; "  e  " ; TAB (30) ; " 

po  " ; T AB ( 45 ) ; "  th  ";TAB<62>;"  beta  ":PRINT  #1, 

3220  FOR  1=1  TO  NPL: PRINT  TAB (9) ; I ; TAB ( 16)  s  s PRINT  USING 
"###.###" :E< I) ; : PRINT  TAB <31 );: PRINT  USING 
"##.##";PO(I)  ;  :  PRINT  TAB  (42)  s:  PRINT  USING 
"###.###"sT(I) ; 

3230  IF  I =NPL  THEN  GOTO  3240  ELSE  PRINT  TAB (60) PRINT  USING 
"###. ###" . BETA ( I ) 

3240  PRINT  #1,  TAB  (9)  ;  I;  TAB  (16)  ;  .-PRINT  #1,  USING 

"###.###••;  E(I)  5  sF'RINT  #1,  TAB  (31 )  5  SPRINT  #1,  USING 
"##.##":PO(I) ; :PRINT  #1,  TAB(42) ; :PRINT  #1,  USING 
”###. ###" ; T ( I ) 5 

3250  IF  I =NPL  THEN  GOTO  3260  ELSE  PRINT  #1,  TAB (60) ? s PRINT 
#1,  USING  " ### . ### " s  BET A ( I ) 

3260  NEXT  I 

3270  PRINTsPRINT  TAB < 5) ; " joi nt  no. " ; TAB ( 1 8) ; "  r  " | TAB (30) j " 
res  " ; TAB (45) ; "  pe  "; TAB (62) 5"  u  "sPRINT 
3280  PRINT  #1, sPRINT  #1,  TAB (5) ; " joi nt  no. " ; TAB ( 18) ; ”  r 

" ; TAB (30) ; "  res  "; TAB (45) 5"  pe  ";TAB(62)s"  u  "1  PRINT 
#1, 

3290  FOR  1=0  TO  NPL-1 sPRINT  TAB (9) 5  I ; TAB ( 16) PRINT  USING 
" ### . ### " { R ( I ) ; sPRINT  TAB (31 ) s IRES(I) ;TAB(42) ; sPRINT 
USING  "###. ###"sPE( I) ; sPRINT  TAB (60) 5 : PRINT  USING 
"###. ###" ; U ( I > 

3300  PRINT  #1,  TAB (9) ; I; TAB (16) 5 sPRINT  #1,  USING 

"###.###" ;R( I) s sPRINT  #1,  TAB (31 ) 5  IRES ( I ) 5  TAB (42) 5 s 
PRINT  #1,  U  SING  "###.###";PE(I> ; sPRINT  #1,TAB<60>? 
sPRINT  #1,  USING  "###. ###" ; U ( I ) 

3310  NEXT  I 

3320  REM  INPUT  "  1  to  cont i nue " ; YCON 
3330  PRINT: PRINT: PRINT: PRINT 

3340  PRINT  #l,s  PRINT  #1, sPRINT  #1, sPRINT  #1,:RETURN 
3350  REM  Soul  t  ion  o-f  Annular  buckling  Equation  .The  D.E. 
being  solved  is: 


REM  D.  E.  =ri  "'"4*  (d4f  )  +2*r  i  •""3*  (d3T  )  -ri"'" 2*  ( l+2*nc"  2+ 
betac(i)*(alfa(i)  '2)  -r  i  "'•2*al  -f  a  ( i  )  •'"2)  *  (d2f  )  + 
r  i  *  ( l+2*nc  "  2+betac  ( i  )  *al  f  a  ( i  >  "‘"2-r  i  "'•2*al  f  a  ( i  )  '•2)  *  < 
d  1  f  )  -  (nc  '"2*  ( 4-nc  '"'2-*-betac  ( i  )  *al  f  a  ( i  )  ""2 
+ri  '2*al  -f  a  ( i  )  ''2)  *  (dOf  )  =0 
‘  VAR:  RI=r(i-l)  or  r(i) 
di=d'i/dr  'i 

■f=f(r),  NC=no.  o-f  waves  in  circum-ferntial 
direction 

al  -f  a  ( i  )  =sqr  < -so-f  ( i  ) /d  ( i  )  )  i -f  sof(i)<0  (tension) 


3390  '  alf  a  (i  )  =sqr  (sof  (i  )  /d  (i  )  >  if  sof(i)>0 

(compr essi on ) 

d  (i  )=e(i  >*t  (1  )  ••3/12/  (l-po(i  >  '"2) 

betac ( i ) =socf ( i ) /sof  ( i >  betac(i>>0  means 

compressi on 

3400  betac ( i ) <0. means  tension 

betac  (  i  )  *al  f  a  (  i  )  "-2=socf  (i  )  /d  (i  )  =sof  cd  <  i  ) 

3410  '  D.  E.  =ri  ""4*  (d4f  )  +2*ri'"3*  (d3f  )  -ri'"-2*  ( 1 +2*nc"''2+socf  d  <i  )  - 
r  i  ‘"-2*sof  d  (  i  ) ""  2>*<d2f)+ri*(  1 +2*nc'"'2+socf  d  (  i  )  - 
r i  2*sof d ( i ) *  (d If )  -  (nc’"'2*  ( 4-nc  '"’2+socf  d  ( i  )  + 
r i  2*sof  d  <  i  )  )  *  < dOf  >  =0 

3420  •  IF  S0FD(I>.=0  THEN  ALFA < I > =SQR < S0FD  ( I ) ) 

3430  '  Assuming  a  power  series  solutions 

f ( i  , k , r l ) -sum C  j  =0 , nn 1 Ca  < i , k , j ) *r i <  j  +sr  ( i , k ) +i *si  ( i , k > ) 

; C  k  =  l  to  4] 

3440  '  calculation  of  the  roots.  The  indicial  eq.  iss 
I .  E.  =  (s-1 )  •'■'4-  ( 2+2*nc  ''"2+socf  d  >  *  (s-1 >  ■'-2+  (  <nc-l )  ''2+ 
socf d* ( l-nc‘2) ) =0 

3450  ‘  The  roots  are  calculated  in  the  follwing  subroutine 
3460  GAMMA  ( I )  =  1 +NONC+S0CFD  (  I )  /2s 

ALAMDA ( I ) -4*NC*NC+2*NC*NC*S0CFD ( I  > +S0C  FD ( I ) *S0CFD ( I ) 

/ 4 :  REM  al  amda  ( l  )  =  ( 2*NC  " 2+SQCFD  ( I  >  /2)  •••2-4*NC"-2*  ( NC"-2- 1 ) 
3470  IF  ALAMDA ( I >  >0  THEN  GOTO  3480  ELSE  IF  ALAMDA ( I ) =0  THEN 
GOTO  3590  ELSE  GOTO  3640 
3480  ’  <<<<  lamc!a>0  >>>>>> 

3490  PHI  1 «GAMMA < I ) +SQR ( ALAMDA (I) ) s  PHI2=GAMMA ( I ) - 
SOR (ALAMDA ( I ) ) 

3500  IF  ABS  <PHI 1— CINT (PHI  1 ) ) <  =  1E— 10  THEN  PHI 1=CINT (PHI  1 ) 

3510  IF  ABS (PHI2-CINT (PHI 2) ) <«1E— 10  THEN  PHI2=CINT (PHI2) 

3520  IF  PH 1 1>0  THEN  PRE 1 =SQR ( PH 1 1 ) s  P I M 1 =0 

3530  IF  PHI  1=0  THEN  PRE 1=0: PIM 1=0 

3540  IF  PHI  1/0  THEN  PRE 1=0 ! : PIM1=SQR (-PHI  1 ) 

3550  IF  PH I 2,0  THEN  PRE2-SQR (PHI2) s PIM2=0 

3560  IF  PH 1 2=0  THEN  PRE2=0s PIM2=0 

3570  IF  PHI 2< 0  THEN  F'RE2=0  !  s  PIM2=SQR  ( -PHI 2) 

3580  GOTO  3670 

359 0  '  < <  <  <  <  1  a  m d  a = 0  > >  > >  > > > 

3600  PH 1 1 =GAMMA ( I ) s  PH  1 2=GAMMA ( I > 

3610  IF  PHI1>0  THEN  PRE1=SQR (PHI  1 ) s PIM 1=0; 

PRE2=PRE1 s  PIM2=PIM1 

3615  IF  PHI  1=0  THEN  PRE1=0: PIMl=Os PRE2=PRE1 s PIM2-PIM1 
3620  IF  PHI l<0  THEN  PRE 1=0 ! : PIM1=SQR (- 
PH  I  1 ) s  PRE2=PRE 1 : P I M2=P I M 1 
3630  GOTO  3670 

3640  '  < < < < < <  <  1  a  md  a  <  0  >  >  > > > >  >  >  > 

3650  PRE 1=SQR (  (GAMMA ( I ) +SQR (GAMMA ( I >  *GAMMA ( I )  + 

ABS (ALAMDA (I))))/2)s  PIM1=  SQR ( -ALAMDA ( I )> /2/PRE  1 
3660  PRE 2= -PRE 1 s  P I M2=P IM 1 
3670  REM  Derivation  of  the  roots 
3680  SR ( I , 1 ) = 1 +PRE 1 s  S I  < I , 1 ) =PIM1 
3690  SR ( I , 2 ) = 1 -PRE 1 : S I (1,2) =-PIMl 
3700  SR  < I , 3) =l+PRE2s  SI (1,3) =PIM2 
3710  SR (1,4) =l-PRE2s  SI (1,4) =-PIM2 


3720  IF  <PRE2=-F'RE1  AND  PIM2-PIM1)  AND  PIM2O0  THEN 

SR ( I , 1 ) -1+PRE1 i SI < I , 1 ) -PIM1 ! SR < I ,2) -1+PRE1 : SI < I , 2) — 

P I M 1 : SR ( 1 , 3 )  =  1  +PRE2: SI ( I , 3) -PIM2: SR ( I , 4 ) - 1+PRE2: 

SI ( 1,4) =-PIM2 

3730  REM  FOR  J  =  1  TO  4: BIGA-SR ( I  ,  J  )  : FOR  K=J  TO  4 
3740  REM  IF  BIGAOSR  ( I  ,K>  THEN  GOTO  2343  ELSE 

BIGA=SR ( I , K) : SR ( I  ,  K ) =SR ( I , J) : SI < I , K) -SI < I , J ) 

3750  REM  K: SR ( I , J) -BIGAs  NEXT  J 

3760  PRINT: TAB < 10) ; "  i  " j TAB (30) ; "  sr  •• ;  TAB  (45)  ;  "  si 
" : PR  I  NT 

3770  FOR  K= 1  TO  4: PRINT; TAB ( 12) ; I : TAB (30) ; SR ( I , K) ? 

TAB (45) ; SI (I, K) : NEXT  K: PRINT 
3780  *  Calculation  of  the  various  solutions 

3790  K=1  s  IF  SR  <1,1)0  0  OR  SR  <  1 ,  1 )  >0  AND  SI(I,1>=0  THEN  GOTO 
3800  ELSE  GOTO  3810 

3800  GOSUB  6000 s  J 00  < 1 ) - JO : J 02= J 0 / 2 : A  < I , 1 , J 02 ) = 1 : GOSUB 
4650: GOTO  3820 

3810  A ( I , 1 ,0) =1 ! : A < I ,2,0) =0 ! : GOSUB  4260:G0T0  3880 
3820  K=2:  '  -Mi, 2) — will  be  consider  for  s=sr<i,2)  only 

3830  IF  SR < I ,2) < >SR < I , 1 )  THEN  GOTO  3840  ELSE  GOTO  3860 
3840  GOSUB  6000s JOO <2) -JO: IF  J0>0  AND  INDX>0  THEN  GOTO  3870 
ELSE  GOTO  3850 

3850  JO-JOO  < INDX ) : J02-JO/2: A < I , 2 , J02) = 1 : GOSUB  4650s  IF 

ABS  <  F  < I , 2 ) -F ( I , 1 ) ) < 1 E- 1 0  THEN  GOTO  3860  ELSE  GOTO  3880 
3860  JO- JOO < 1 ) /2: A  < I , 2 , JO) =1 s IDS=1 s  KKK=1 : GOSUB  5220s  GOT  3 
3880 

3870  IF  INDX >1  THEN  GOTO  3850  ELSE 

IDS=1 : KKK=INDXs  JO-JOO (KKK) : GOSUB  4840s  GOTO  3880 
3880  K=3s IF  ABS <SI < I ,3) ) >0  THEN  GOTO  4010 

3890  IF  SR (1,3) =SR  < I , 1 )  AND  SR < I , 3) =SR < 1 , 2)  THEN  GOTO  3990 
ELSE  IF  SR (1,3) =SR (1,1)  THEN  KKK— 1  ELSE  IF 
SR ( I ,3) =SR ( I ,2)  THEN  KKK— 2  ELSE  GOTO  3910 
3900  GOTO  4000 

3910  GOSUB  6000s JOO (3) -JO: IF  J0>0  AND  INDX>0  THEN  GOTO  3950 
ELSE  GOTO  3920 

3920  JO- JOO <  KKK )/2sA(I,3,0)*l!s GOSUB  4650 : I F  ABS  <  F  <  I  ,  3  >  - 
F ( I , 1 ) ) <  =  1 E- 1 0  THEN  KKK- 1 s  GOTO  3980 
3930  IF  ABS ( F ( I , 3 ) -F ( I , 2 ) ) <  =  1 E— 1 0  THEN  KKK=2s  GOTO  3980 
3940  GOTO  4040 

3950  IF  IND ( 1 ) >0  THEN  KKK- 1: GOTO  3980 
3960  IF  I ND ( 2 ) >0  THEN  KKK=3:GCT0  3980 
3970  IF  IND (3) >0  OR  IND(4)>0  THEN  KKK=3:G0T0  3920 
3980  JO-JOO (KKK) /2: A(I,3,J0)=1: IDS-1 s GOSUB  4840s  GOTO  4040 
3990  KKK-1 : JO-JOO  < 1 ) /2s  A(I,3,0)  =  1: IDS-2: GOSUB  5220: GOTO  4040 
4000  JO-JOO (KKK) /2s A< I ,3, JO) =1 s IDS-1 : GOSUB  5220s  GOTO  4040 
4010  IF  SR (1,3) -SR (1,1)  AND  SI ( I ,3) -SI < I , 1 )  THEN  GOTO  4030 
4020  A< I ,3,0) -Is A(I , 4 , 0) =0s GOSUB  4260: GOTO  4220 
4030  A ( 1 , 3 , 0  >  = 1 s  A ( 1 , 4 , 0 ) =0  s  GOSUB  4260  s  GOTO  4220 
4040  K=4:  'f(i,4  will  be  considered  only  for  s=sr(i,4)  only 

4050  IF  (SR  < I , 4) <  >SR (1,1))  AND  (SR ( I ,4) < >SR < I ,2) )  AND 
( SR  < I , 4  >  <  >SR (1,3) )  THEN  GOTO  4110 
4060  IF  SR (I, 4) -SR (1,1)  AND  SR < I , 4) -SR ( I , 2)  AND 
SR (1,4) -SR (1,3)  THEN  GOTO  4180 


4070  IF  SR (1,4) =SR  < 1 , 1 )  THEN  KKK-1  ELSE  IF  SR ( I , 4 > =SR ( I , 2) 

THEN  KKK=2  ELSE  IF  SR  <  I  , 4)  =>SR  <  1 , 3)  THEN  KKK-3  ELSE  SOTO 
4090 

4080  SOTO  4200 

4090  IF  <SR<I,4)=SR(I,1>  AND  SR ( I , 4 ) =SR ( I , 2) )  THEN  KKK=1 

ELSE  IF  ( SR (1,4) =SR (1,3)  AND  SR ( I , 4 ) =SR ( I  ,  1 ) )  THEN  KKK= 1 
ELBE  IF  ( SR (1,4) =SR (1,2)  AND  SR ( I , 4 ) =SR ( I , 3 ) >  THEN  KKK-2 
ELSE  GOTO  4220 
4300  GOTO  4190 

4110  GOSUB  6000: J00<4) =J0: IF  J0>0  AND  INDX>0  THEN  GOTO  4120 
ELSE  4160 

4 1 20  FOR  1 1 1  —  1  TO  3 : . F  I ND  < 1 1 1 >  >0  THEN  KKK= III: GOTO  4 1 50 
4130  NEXT  I  I  I 

4140  IF  I ND  <  4 )  >0  THEN  KKK=4: GOTO  4 1 60 
4150  A ( I , 4 , 0 ) = 1 ! : I DS= 1 : GOSUB  4840: GOTO  4220 

4160  J0«J00 (4) /2: A ( I ,4, JO) =1 ! : GOSUB  4650:F0R  111  =  1  TO  3:  IF 
ABS  <F  ( I  ,4)  -F  ( I  ,  1 1 1 )  )  O1E-10  THEN  KKK* III:  GOTO  4210 
4170  NEXT  1 1 1: GOTO  4220 

4190  J0=J00 ( 1 ) /2: A (I ,4, JO) =1 ! : IDS=3: KKK=1 : GOSUB  5220: GOTO 
4220 

4190  JO=JOO (KKK) /2: A ( I , 4, JO) — 1 ! : IDS=2: GOSUB  5220: GOTO  4220 
4200  JO=JOO (KKK) /2: A< I ,4, JO) “1 ! : IDS=1 : GOSUB  5220: GOTO  4220 
4210  J  0= J  00  ( KKK )  /  2 :  A  ( I  ,  4 ,  J  0 )  =  1  1  :  I  DS=  1 :  GOSUB  5220.-GCT0  4220 
4220  IF  YEPR$="no"  THEN  RETURN  ELSE  PRINT; TAB (5) ; "  l 

" j  TAB ( 20 ) ; ”  fi  ";TAB<35>;"  dfi  " ; TAB (50) ; "  d2  i 
"sTABCSS);"  d3fi  ": PRINT 

4230  FOR  K-l  TO  4: PRINT; TAB (5) ; I » TAB ( 17) ; F ( I , K> ; TAB (30) ; 

DF ( I ,K) ; TAB (47) ; D2F ( I  ,K);T  AB (65) ; D3F ( I , K) : NEXT 
K: PRINT 

4240  REM  FOR  K=1  TO  4:F'RINT  #2 ,  "  i  ,  k  ,  F  ( I  ,  K)  ,  D2F  ( I  ,  K )  , 

D3F ( I , K ) ,d4f ( i , k ) 

I , K , F  < I , K )  , DF  < I , K  > ,D2F(I,K) ,D3F ( I , K ) , D4F ( I , K ) : 

NEXT  K: PRINT 
4250  RETURN 

4260  REM  Solution  with  s=sr+i*si (Reccurrence  formula) 

4270  NNS2=2*NNS: FOR  J=2  TO  NNS2  STEP  2:J2=J/2 
4280  ANOMR 1 =J-2+SR ( I , K ) : ANOM 1 1 =S I  < I , K ) 

4290  ANOMR=ANOMR 1 *ANOMR 1 -ANOM 1 1  * ANOM 1 1 - 
NONC :  ANOM  I  =2*AN0MR  1  *  ANOM  1 1 
4300  IF  <AN0MR=0  AND  AN0MI*0)  OR  SOFD ( I ) *0  THEN 
A ( I , K , J 2 ) =0 : A ( I , K+ 1 ,J2)=0:G0T0  4380 
4310  DEN0MR2=J-1+SR ( I ,K) :DEN0MI2-SI (I,K> 

4320  T I R 1 ®DENQMR2*DEN0MR2- 

DENOM I2*DEN0M 12: Till =2*DEN0MR2*DEN0M 1 2 
4330  TIR2»TIR1*TIR1-TIU*TII1:TII2=2*TI  1 1  *T I R 1 
4340  DENOMR=T I R2-2*GAMMA ( I >  *TIR1+ ( 1-NC*NC) *  (  1- 

NC*NC+S0CFD( I > ) : DENOMI=TI  I2-2*GAMMA ( I ) *TI I  1  : 
DNOR=DENOMR*DENOMR+DENOMI#DENOMI 
4350  ARR=-SOFD ( I ) *  <DENOMR*ANOMR+DENOMI*ANOMI ) /DNOR 
4360  All =-SOFD ( I ) * ( DENOMR*ANQM I -DENOM I *ANOMR ) /DNOR 
4370  A ( I , K , J 2 ) =ARR*A ( I , K , J2- 1 ) -A 1 1 *A ( I , K+ 1,J2- 

1) :A(I ,K+1,J2)=AII*A<I ,K,  J2-1 ) +ARR*A ( I , K+l , J2-1 ) 

4380  NEXT  J 

4390  K=K: JMAX=NNS:KPP=K: J0=0: GOSUB  5370 


4400  SU  <  1  > “SUM: DSU  < 1 ) =DSUM:  D2SU ( 1 ) “D2SUM: D3SU ( 1 >  =D3SUM: 
D4SU  < 1 ) “D4S  UM 

4410  K=K+1 :  JMAX“NNS;  KPP“K:  J0=0:  GOSUB  5370 
4420  SU ( 2 ) =SUM : DSU  <  2 ) “DSUM : D2SU  <  2 ) “D2SUM : D3SU ( 2  > “D3SUM : 
D4SU (2) “D4S  UM 

4430  K“K- 1 : G ( 1 > “COS  < S I < I ,K> *LOG (RI > > : 

(3  (2)  “SIN  (SI  (I,K)*LOB(RI)  > 

4440  DG ( 1 ) =-S I (I,K)*G(2> /RI:DG(2)=SI < I ,K> *G < 1 > /RI 
4450  D2G ( 1 ) =-SI ( I , K ) *DG(2) /RI+SI ( I ,K> *G (2) /RI/RI 
4460  D2G (2) “SI (I ,K)*DG< 1 ) /RI-SI ( I ,K) *G < 1 ) /RI/RI 
4470  D3G ( 1 ) =-SI <I,K>*D2G<2> /RI+2*SI <I,K)*DG(2)/RI/RI- 
2*SI(I,K)*G(2)/RI/  RI/RI 
4480  D3G (2) “SI ( I ,K) *D2G ( 1 ) /RI- 

2*S1 (I ,K)*DG<1) /RI/RI+2*SI ( I ,K> *G ( 1 > /RI/R  I/RI 
4490  D4G ( 1 ) “-SI ( I , K) *D3G (2) /RI+3*SI ( I , K) *D2G (2) /RI/RI- 

6*SI (I ,K)*DG<2> /R  I/RI/RI+6*SI (I ,K) *0(2) /RI/RI /RI/RI 
4500  D4G (2) “SI < I ,K) *D3G ( 1 ) /RI- 

3*SI <I,K)*D2G(1>/RI/RI+6*SI ( I ,K> *DG ( 1 > /RI  /RI/RI- 
6*SI < I, K)*G<1) /RI/RI /RI/RI 
4510  HS=Oi FOR  IH=1  TO  2: FOR  JJH=1  TO  2:HS“HS+1 
4520  H (HS) =SU ( IH) *G ( J JH) 

4530  DH (HS) =DSU ( IH) *G  (JJH) +SU(IH>  *DG (JJH) 

4540  D2H (HS) =D2SU ( IH) *G ( J JH) +2*DSU ( IH) *DG ( JJH) + 

SU ( IH) *D2G (JJH) 

4550  D3H  <HS) “D3SU ( IH) *G (JJH) +3*D2SU ( IH) *DG ( JJH)  + 

3*DSU ( IH) *D2G (JJH)  +SU ( I  H)*D3G(JJH) 

4560  D4H (HS) “D4SU ( IH) #G ( JJH) +4*D3SU ( IH) *DG (JJH) + 

6*D2SU ( IH) *D2G (JJH  ) +4*DSU ( IH) *D3G ( J JH) + 

SU( IH)*D4G(JJH) 

4570  NEXT  JJH: NEXT  IH 

4580  F ( I , K+l ) “H ( 1 ) -H (4)  :  F ( I ,K> “H (2) +H (3) 

4590  DF  < I , K+ 1 ) =DH ( 1 ) -DH ( 4 ) : DF ( I , K ) =DH ( 2 ) +DH ( 3 ) 

4600  D2F ( I , K+ 1 ) =D2H ( 1 ) -D2H ( 4 ) : D2F ( I , K > “D2H ( 2  >  +D2H ( 3  > 

46 1 0  D3F ( I , K+ 1 ) =D3H ( 1 ) -D3H ( 4 ) : D3F ( I , K ) “D3H ( 2 ) +D3H ( 3 ) 

4620  D4F ( I , K+l ) “D4H < 1 ) -D4H ( 4 ) : D4F ( I , K ) =D4H ( 2 ) +D4H ( 3 ) 

4630  GOSUB  6290: K-K+l : GOSUB  6290: K=K- 1 : IF  YEPRf *"no"  THEN 
RETURN  ELSE  PRINT" k , k+l , f ( i , k >, i ( i , k+l ) 

";K,K+1 ,F(I,K) , F  < I , K+l > 

4640  RETURN 

4650  REM  Reccurrence  formula  and  solution  for  s=sr 
4660  NNS2“2*NNS: FOR  J=2+J0  TO  NNS2  STEP  2:J2«J/2 
4670  AN0M1=J+SR(I  ,K)-2:  ANOM=ANDM1*ANOM1-NONC:  IF  AM0M=0  OR 
SOFD ( I ) “0  THEN  A ( I , K , J2 > “0: GOTO  4710 
4680  SS“J+SR ( I ,K) : SS1=J+SR  < I ,K> -  1 : SS2=J+SR ( I , K ) - 
2: SS3=J  +SR ( I ,K) -3: SNN= 1 +2*NC*NC+S0CFD < I ) : 
DEN0M“SS*SS1*SS2*SS3+2*SS*SS1*SS2-SNN*SS*SS  1+SNN*SS- 
NC*NC*  (4-NONC+SOCFD  <  I )  > 

4690  IF  A ( I , K , J2- 1 ) <  >0  THEN  A(I,K,J2)=- 

SOFD ( I ) *ANOM/DENOM*A (I , K , J2-1 >  ELSE  A(I,K,J2)*0 
4700  REM  DEN0M1“J+SR(I ,K)-1:DEN0M=<DEN0M1*DENCM1- 

PHI 1 )  •*  (DEN0M1*DEN0M1-PHI2>  :  IF  A  ( I  ,  K ,  J2-1 )  <  >0  THEN 
A ( I , K , J2> “-SOFD ( I ) *aNOM/DENOM*A ( I , K , J2- 1 )  ELSE 
A ( I  , K , J  2  > “0 


4710  NEXT  J : REM  PRINT  ANOM , DENOM* PRINT 

"  i  =  "  ;  I ,  "  k  =  "  ;  K ,  "  j  =  "  s  J  ,  "  a  ( i  ,  k. ,  j 2 >  “ "  j  A  ( I ,  K ,  J 2 )  :  NEXT  J 
4720  GOTO  4830s  REM  PRINT 

#2, "a ( " s I ; " , " ; K* " , " j  JO" ) =" *  A ( I , K , JO) s GOTO  4080 
4730  JMAX*NNSs GOSUB  5370 1  RETURN 

4740  NNS2=2*NNSs  SUM=0«  DSUM=0t D2SUM=0s D3SUM=0:  D4SUM=0s FOR  J=J0 
TO  NNS2  STEP  2sJ2=J/2 
4  750  SUM=SUM+A  ( I  ,K, J2>  *RI  "' ( J+SR  <  I  ,KPP>  ) 

4760  DSUM=DSUM+A(  I  ,K,  J2)  *  (J+SR  <  I  ,KPP>  )  *RI""  (J+SR  ( I  ,KPP)  -1 ) 

4770  D2SUM=D2SUM+A ( I , K , J2 ) * ( J+SR ( I , KPP ) ) *  <  J+SR ( I , KPP ) - 
1)*RI "(J+SR (I ,KPP  ) -2 > 

4780  D3SUM=D3SUM+A  ( I ,  K ,  J  2  >  *  (J+SR  ( I ,  KPP  >  )  *  ( J  +SR  ( I  v  KPP  >  - 
1 )* (J+SR ( I , KPP) -2  ) *RIA (J+SR( I , KPP) —3) 

4790  D4SUM=D4SUM+A ( I , K , J2) * ( J +SR ( I , KPP ) ) * (J+SR ( I , KPP ) - 

1 )  *  (J+SR  ( I  ,KPP)  -2  )  *  (J+SR  ( I  ,  KPP )  -3 )  *R  I  (J+SR  ( I , KPP)  -4) 
4800  REM  PRINT  " k , kpp , SUM , DSUM , D2SUM , D3SUM , J " ; K ,KPP , SUM , 
DSUM , D2SUM , D3SUM , D4S  UM,J 
4810  NEXT  J 
4820  RETURN 
4830  KPP*Ks GOSUB 

4 730 : F ( I , K ) =SUM : DF  < I , K >  =DSUM s  D2F ( I , K > =D2SUM : D3F ( I , K ) =D3SUM : D4 
F ( I , K  ) =D4SUM: GOSUB  6290s  IF  YEPR*="no"  THEN  RETURN  ELSE 
PRINT  " F ( I , K ) , K  " 5  F ( I , K ) , K s  RETURN 
4840  REM  Reccurrence  -formula  and  solution  using 
D-  ndsF/DSAnds.  Nds-no.  of  differentiation 
4850  NNS2ss2*NNSs  FOR  J=»2+J0  TO  NNS2  STEP  2:J2=J/2 
4860  AN0M1=J+SR(I ,KKK)-2s AN0M»ANQM1*AN0M1-NC*NC: IF  AN0M=0  OR 
SOFD ( I ) *0  THEN 

AAA ( J2) =0s  DA ( J2>  =0: D2A ( J2) =0s  D3A ( J2) =0: GOTO  4930 
4870  DN0M=2*N0M 1 s  D2N0M=2s  D3N0M=0 
4880  DEN0M2* J  +SR  ( I  ,  KKK)  -1 :  DEN0Ml=DEN0M2'  -2- 
GAMMA ( I ) s  DENOM= 1 / <  DENOM 1 ' 2- AL  AMDA ( I ) ) 

4890  DDEN0M=2*DEN0M  1  *2*DEN0M2  s  D2DEN0M=4*  ( 2-»DEN0M2  -2+DEN0M  1 ) 
s  D3DEN0  M=4*<  4*DEN0M2+2*DEN0M2) s AAA ( J2> =NCM*DENOM 
4900  DA(J2)  =  (DN0M*DEN0M-N0M*DDEN0M*DEN0M"-2)  s  IF  IDS-=1  THEN 
GOTO  4930 

4910  D2A(J2)=(D2N0M*DEN0M-2*DN0M*DDENQM*DEN0M'2- 

NOM*D2DENOM*DENOM""’2+2*N  OM*DDENOM-"'2*DENOM"",3 >  s  IF  IDS«2 
THEN  GOTO  4930 

4920  D3A  ( J2>  =  (D3N0M*DEN0M-3*D2N0M*DDEN0M*DEN0M  '"2- 
3*DN0M*D2DEN0M*DEN0M-'-2 

+6*DN0M*DDEN0M"-2*DEN0M-  3+6*N0M*D2DENQM*DDEN0M*DEN0M""-3- 
N0M*D3D  ENGM*DEN0M'-2-6*N0M*DDEN0M--3*DEN0M '4 ) 

4930  NEXT  J s REM  PRINT  " j , nom ,noml , denom , da ( j2> , aaa ( j2) 

"5 J,N0M,N0M1 , DENOM, DA <J2> ,AAA(J2) sNEXT  J 
4940  J02=J0+2s J022=J02/2s A < I , K , J022) =-SOFD ( I ) *DA < J022) s IF 

IDS=2  THEN  A ( I , K , J022 ) «-SOFD ( I ) *D2A ( J022)  ELSE  IF  IDS=3 
THEN  A ( I » K , J 022 ) D3 A  <  J  022 ) *SOFD ( I ) 

4950  NNS2=2*NNSs  FOR  J=4+J0  TO  NNS2  STEP 

2:  J2ssJ/2s  SUMD=Os  SUM2D=0s  SUM3D==0s  FOR  LL=J0+2  TO  J  STEP 
2sLL2=LL/2s IF  AAA(LL2>=0  AND  DA(LL2>=0  AND  D2A(LL2)=0 
AND  D3A (LL2) =0  THEN  GOTO  5030 
4960  SUMD=SUMD+DA (LL2) /AAA (LL2) s IF  IDS=1  THEN  GOTO  4990 
4970  SUM2D**SUM2D+D2A  (LL2)  /DA  (LL2)  s  IF  IDS=2  THEN  GOTO  4990 


4980  SUM3D=SUM3D+D3A(LL2> /D2A (LL2) 

4990  NEXT  LL 

5000  A(I,K,J2)=SUMD*A(I,KKK,J2> s IF  IDS-1  THEN  B0T0  5030s' 
a=da/ds 

5010  A  < I , K , J2> =SUM2D*A  < I , KKK , J2>  s IF  IDS=2  THEN  GOTO  5030s' 
a=d2a/d2s 

5020  A  ( I  ,  K ,  J  2 ) =SUM3D#A < I , KKK , J2)  s  '  A 

=D2A/D2S 

5030  NEXT  JsREM  PRINT 

"da  <  j2)  ,  aaa  ( j2)  , j";DA(J2)  ,AAA(J2)  ,JsPRINT 
"  i  =  "  5  I  ,  "  k  =  "  ;  K  ,  "  j  =  "  ;  J  ,  "  a  ( i  ,  k  ,  j  2  >  =  "  ;  A  <  I  ,  K  ,  J  2 ) sNEXT  J 
5040  J 02= JO/2 s  A ( I , K , JO) =0 s  KPP=KKK s  GOSUB  4730 
5050  IF  RIO.  000001  THEN  B0=0  ELSE  B0=LCG(RI)  'IDS 
5060  IF  ROC.  000001  OR  IDS-1  <0  OR  LQG(RI)=0  THEN  BOO  ELSE 
BO-LOG  (RI)  "-(  IDS-1) 

5070  IF  ROC.  000001  OR  IDS-2C0  OR  L0G(RI)=0  THEN  B2=0  ELSE 
B2=LQG<RI  )-'•(  I DS-2) 

5080  IF  RIO.  000001  OR  IDS-3C0  OR  L0G(RI)=0  THEN  B3=0  ELSE 
B3=L0G(RI  )•-••(  I DS-3) 

5090  IF  RIO.  00000 1  OR  IDS-4C0  OR  L0G<RI)=0  THEN  B4=0  ELSE 
B4=L0G  ( R I )  •"•  ( I DS-4 ) 

5 1 00  F ( I , K ) =F  < I , KKK ) *BO+SUM 

5110  DF ( I , K ) ~DF < I , KKK ) *B0+DSUM+F ( I , KKK > * I DS*B 1 / R I 
5120  D2G=IDS*  ( IDS-1 )  +B2/RI  /RI+IDS*B1*  (-1/RD/RI 
5 1 30  D3G= I DS*  < I DS- 1 ) * ( I DS-2 ) *B3/R I /R I /R I -3* I DS* ( I DS- 
1)*B2/RI/RI/RI+IDS*  B1*2/RI''3 
5140  D4G=IDS*  (IDS-1 )  *  ( IDS-2)  *  <  IDS-3)  *B4/RI '-4-6*IDS*  ( IDS- 
1 )  *  <  IDS-2)  *B3/R  I  '4+1 1*IDS*  (IDS-1 )  *B2/RI''4- 
IDS*B1*6/RK4 

5 1 50  D2F  ( I  ,  K )  =D2F  ( I  ,  KKK )  *L0G  (RI)-'I  DS+2*DF  ( I  ,  KKK ) 
*IDS*B1/RI+F(I,KKK  )  *D2G  +D2SUM 
5 1 60  D3F ( I , K ) =D3F ( I , KKK ) *L0G ( R I ) - I DS+3+D2F ( I , KKK ) 

*IDS*B1 /R1+3+DF ( I  , KKK )  +D2G+F ( I , KKK) *D3G+D3SUM 
5 1 70  D4F ( I , K ) =D4F ( I , KKK >  *L0G  ( R I )  - 1 DS+ 

4+D3F ( I , KKK) * IDS+B 1 /R I +6+D2F  (  I , KKK 
) *D2G+4*DF ( I , KKK ) +D3G+F ( I , KKK  >  +D4G+D4SUM 
5180  GOSUB  6290 

5190  IF  YEPR*="no"  THEN  RETURN 

5200  PRINT  " K , F  < I , K )  , DF ( I , k ) , D2F ( I , K ) , D3F ( I , K ) , 

d4-f  ( l  ,  k  )  ,  KKK  ,  BO ,  B 1  ,  B2 ,  B3 ,  B4  ,  DG  ,  D2G ,  D3G,  d4g  "  s  K 
,F(I,K) , DF ( I , K ) , D2F ( I , K ) , D3F ( I ,K) ,D4F(I,K) , KKK 

,B0,B1 ,B2,B3,B4, DG , D2G , D3G , D4Gs PRINT  " - RI=";RI 

5210  RETURN 

5220  REM  Solution  for  equal  roots 

5230  NNS2=2*NNS: FOR  J-0  TO  NNS2  STEP  2s  J2=J/2 

5240  DAA 1 =A ( I , KKK , J2) + ( J+SR ( I , KKK ) - 

1 )  * ( 4* ( J  +SR ( I , KKK ) )  +  ( J  +SR ( I , KKK ) -2 ) -4*NC*NC-S0CFD ( I ) *2 ) 
5250  IF  J=0  THEN  GOTO  5260  ELSE  DAA2=2*S0FD(I)*(J+SR(I,KKK)- 

2)  *A  < I , KKK , J2- 1 >  s  DA  1 =DAA 1 +DAA2 

5260  ANOM 1 = J  +SR ( I , K ) -2 : AN0M=AN0M 1 *AN0M 1 -NC*NC 
5270  SS=J  +SR ( I , K) sSSl=J+SR ( I ,K) -1 : SS2-J+SR ( I , K ) - 
2:  SS3=J+SR ( I ,K) -3: SNN= 1 +2*NC*NC+S0CFD ( I ) : 
DEN0M=SS*SS1+SS2*SS3+2*SS*SS1*SS2-SNN*SS*5S  1+SNN*SS- 
NC*NC* (4-NC+NC+SOCFD ( I ) ) 


5280  IF  J >2  THEN  GOTO  5330 

5290  IF  J=0  AND  DAA1=0  THEN  GOTO  5340  ELSE  IF  J=0  AND  DAAIOO 
THEN  PRINT  "TRY  ADDITIONAL  TERMS" s BEEP: BEEP: BEEP: STOP 
53o0  IF  J=2  AND  ( DEN0M=0  AND  AN0M*0>  THEN  PRINT  "try 
addi ti onal  terms" : BEEP: BEEP: BEEP: STOP 
5310  IF  J=2  AND  DEN0M=0  THEN  A(I,K,0>=- 
D A 1 / SOFD ( I ) / ANOM :A(I ,  K ,  J  2 )  “  1 ! 

5320  IF  J*>2  AND  DENOMOO  THEN  A  ( I ,  K ,  0 )  =0 :  A  ( I  , K , J  2 ) =- 
DA1/DEN0M 

5330  A ( I , K , J2) *— (DA1+A ( I , K , J2-1 ) *S0FD ( I ) *ANOM) /DENOM 

5340  NEXT  J 

5350  KPP=K: GOSUB  4730 

5360  GOSUB  5050: RETURN 

5370  REM  Sub.  Evaluate  a  polynom  at  r=ri  using  nested 
procedure 

5380  REM  VAR: jmax , jO, a ( i , k , j2) , k , kpp 

5390  SUM=0: DSUM=0: D2SUM-0: D3SUM=0: D4SUM*0: J J«JMAX 

5400  JJ2=JJ*2: AO=A (I ,K, JJ> : SJS=J J2+SR < I , KPP) s 

DAO*SJS*AO: D2A0=DA0* (  SJS-1  ) : D3A0=D2A0* (SJS- 
2> : D4A0-D3A0* (SJS-3) 

5410  SUM»SUM*RI*RI+AO: DSUM=DSUM*RI*RI+DAO: 

D2SUM=D2SUM*R I #R I +D2 AO I D3SUM=D3SUM*R I *R I +D3A0 : 
D4SUM=D4SUM*R I *R I +D4A0 

5420  REM  PRINT  " k , kpp , SUM , DSUM, D2SUM , D3SUM , d4sum , J j " ; 

K , KPP , SUM , DSUM , D2SUM , D3  SUM,D  4SUM,JJ 
5430  IF  JJ=<0  THEN  GOTO  5450 
5440  J J=J J-l : GOTO  5400 

5450  IF  SR  < I , KPP) <0  THEN  SP0»1 /R I" ABS (SR < I , KPP) >  ELSE 
SPO«R I 'SR ( I , KPP ) 

5460  SUM*SPO*SUM:  DSUM=SPO*DSUM/RI :  D2SUM=SF'0*D2SUM/RI /RI 

:D3SUM=SP0*  D3SUM  /RI/RI/RI: D4SUM=SP0*D4SUM/R I /RI/RI/RI 
5470  RETURN 

5480  REM  data  for  a  run 

5490  'DATA  1,.2,1 

5500  '  DATA  1 , . 2 , 1 

5510  '  DATA  1,.2,1 

5520  '  DATA  1,.2,1 

5530  '  DATA  4,3,2, 1 

5540  '  DATA  2,0,0, 1 ,0,0, 0,0 

5550  '  data  for  two  part  plate' 

5560  DATA  "aaa" 

5570  DATA  2,30,6 

5580  DATA  2,2 

5590  DATA  1 ,. 3000000 , 1 . 

5600  DATA  1 , . 3000000 , 1 . 

5610  DATA  1 . 0000 1 , . 500 

5620  DATA  2 , 1 , 0 , 0 

5630  DATA  0,10 

5640  REM  Invert  of  a  matrix 

5650  'Var.:  aA(n)=Input  matrix  destroyed  in  computation  and 
replaced  by  resultant  invrese  (columnwise) 

5660  REM  N«  order  of  matrix  A. (not  columnwise) 

5670  '  det*  resultant  determinant 

5680  '  L=  WORK  VECTOR  OF  LENGTH  N 


5690  '  M=  WORK  VECTOR  OF  LENGTH  N 

5700  ’  dim  a  (n )  , r  1  (n >  , rm  (r, ) 

5710  DET =  1  !  : NK=-N 

5720  FOR  K=  1  TO  Ns  NK=NK+N:  L  (K)  =Ks  M  <K)  =K:  KK=NK+Ks  BIGA=AA  (KK) 
5730  FOR  J=K  TO  Ns  I Z=N* < J-l ) a  FOR  I=K  TO  NsIJ=IZ+I 
5740  IF  (ABS(BIGA)-ABS(AA(IJ) ) > >=0  THEN  GOTO  5750  ELSE 
BIGA=AA ( I J ) sL(K)=IsM(K)=J 
5750  NEXT  I: NEXT  J 

5760  J=L  (K>  :  IF  (J-KX-0  THEN  GOTO  5780  ELSE  KI=K~N 
5770  FOR  1=1  TO  N: KI=KI+Ns HOLD=-AA <KI > a J I=KI- 
K+J  s  AA ( K I >  =AA ( J I ) : AA ( J I ) =HOLD: NEXT  I 
5780  I=M  (K)  :  IF  (I-K)OO  THEN  GOTO  5800 
5790  JP=N* ( I - 1 ) : FOR  J=1  TO  N: JK=NK+J s  J I =JP+J : HOLD=- 
AA(JK) s AA<JK)=AA(JI) j AA ( J I > =HOLDs NEXT  J 
5800  IF  BIGAOO  THEN  GOTO  5810  ELSE  DET=0 !: RETURN 
5810  FOR  1  =  1  TO  Ns  IF  <I-K>=0  THEN  GOTO  5820  ELSE 
IK=NK+I s AA  < IK) =AA ( IK) / (-BIGA) 

5820  NEXT  I 

5830  FOR  1  =  1  TO  Ns IK=NK+I s  HOLD=AA (IK) s I J=I-N 
5840  FOR  J=1  TO  Ns  I J  =  I J +Ns IF  (I-K>=0  THEN  GOTO  5850  ELSE  IF 
( J-K) =0  THEN  GOTO  5850  ELSE  KJ=IJ- 
I +K  s  A A ( I J  > =HOLD*AA  <  K J ) +AA  < I J ) 

5850  NEXT  JsNEXT  I 

5860  KJ=K-N:  FOR  J  =  1  TO  NaKJ=KJ+NsIF  (J-K)OO  THEN 
A  A  ( K  J )  =  A  A  ( K  J  >  /  B I G  A 
5870  NEXT  J 
5880  DET=DET*BIGA 
5890  AA (KK) =1 /BIGA 
5900  NEXT  K 
5910  K=N 
5920  K=K- 1 

5930  IF  K=< 0  THEN  GOTO  5990  ELSE  I=L(K> 

5940  IF  ( I — K ) <  =0  THEN  GOTO  5960  ELSE  JQ=N* (K-l ) s JR=N* ( 1-1 ) 
5950  FOR  J=1  TO  NsJK=JQ+JsHOLD=AA(JK)sJI=JR+J:AA(JK)=- 
AA ( J I ) s  AA  < JI >  =HOLD: NEXT  J 

5960  J=M (K> s IF  (J-K)<=0  THEN  GOTO  5920  ELSE  KI=K-N 
5970  FOR  1=1  TO  NsKI=KI+NsHOLD=AA(KI)sJI=KI-K>JsAA(KI)=- 
AA  < J I ) s  AA  < J I ) =HOLDs  NEXT  I 
5980  GOTO  5920 
5990  RETURN 

6000  REM  Sub.  for  checking  case  of  denom=0(s=sr  only) 

6010  '  VARs  prel , pre2 ,  h  j  1  , h  j2 , h  j3 ,h  j4 ,  .jh  1  ,  jh2,  jh3 ,  jh4 , 
indl,ind2,ind3,ind4 

6020  HJ  < 1) =1-SR(I ,K)-PRE2s JH ( 1 ) =CINT (HJ < 1 ) ) s HJ  <2>  =  1- 
SR(I ,K)+PRE2s JH(2)=  CINT(HJ (2) > 

6030  HJ (3) =1-SR ( I ,K) -PRE1 s JH (3) =CINT  <HJ (3) ) sHJ(4)=l- 
SR ( I , K ) +PRE 1 s  JH  < 4 ) =  CINT(HJ<4>  > 

6040  FOR  1 1 1  =  1  TO  4s  NQM2=  <2-SR  ( I  ,K)  -2)  "'2-NC",2 
6050  IF  ABS (HJ ( 1 1  I ) -JH (III) ) <  =  1E— 08  AND ( JH ( 1 1 1 > >0  AND 
JH  ( 1 1 1 )  /2*INT  (JH(III)/2)  )  AND  NCM2O0  THEN 
JO=JH ( 1 1 1 ) s  L ( 1 1 1 ) =JH ( 1 1 1 ) s IND  < 1 1 1 >  =  1 1 1 s GOTO  6070 
6060  I ND (1 1 1 ) =0  s  J0=0  s  L ( 1 1 1 ) =0 
6070  NEXT  III 


6080  IF  I ND  (IX  >0  OR  IND<2)<>0  OR  IND  (3)00  OR  IND  (4)  0-0  THEN 
OOTO  6090  ELSE  GOTO  6110 

6090  INDX=1:BIGA=L<1)  :  FOR  1 1 1  =  1  TO  4i  IF  BIGACLUII)  THEN 
BIGA=L (III): INDX=I I 
6100  NEXT  Ills JO=BIGA 
6110  IF  YEF'R#=”no"  THEN  RETURN 

6120  PRINT  JO, INDX , IND  <  1 )  , IND (2)  , IND (3) , IND (4) : RETURN 
6130  REM  calculation  of  internal  forces  at  r=ri 
6140  '  var:  f  ( i , k ) ,df ( i , k ) , d2f ( i , k ) ,d3f ( i , k ) , r i 
6150  DW=DF ( I , K) 

6 1 60  BM=- ( D2F ( I . K > +P0 ( I ) * ( DF ( I , K ) -NC*NC*F  < I , K ) /R I ) /R I >  *D  < I ) 

6 1 70  SH—  <  D3F  (I ,  K )  +D2F  ( I  ,  K )  /  R I  ( NONC+ 1 )  *DF  ( I  ,  K )  /  R I  /  R I  + 
2*NC*NC*F(I ,K> /R  I /RI /RI-INDE*NC*NC* ( 1-PO < I ) ) 
*(DF(I,K)/Ri/RI-F(I,K)/RI/RI/RI)>  *D ( I ) 

6180  IF  YEPR$="no"  THEN  RETURN  ELSE  PRINT 

" i  =  "  ;  I , "k=" j  K , "ri  =  " $  RI s  PRINT" f ( i , k > , dw,bm , sh " ; 

F ( I , K ) ,DW,BM,SH 
6190  RETURN 

6200  REM  E;;planatiom  of  the  signs  for  the  varios  B.C.  ,inc. 

the  needed  equati onl i st  95-130 
6210  REM  ibcw(i)  meaning  graphic  not. 

equations  -  -  - 


6220  REM  0  S.S  E. 

WR=0 , BMR=0 

1  FIXED  E.  ! - 

WR=0,DWR=0 

2  F. FIXED  E.  !  ! - 

WR-0, SHR-0 

6230  REM  3  E.S.S  E.  - 

BMR=0 , SHR-KV ( I ) *WR-0 

4  FREE  E.  - 

BMR-0 , SHR=0 

6240  REM  10  CONTI NOUS  J.  - 

WRL=0 , DWRL-0 , BMRL-0 , SHRL=0 

11  S.S  J.  - ■- - 

WL=0 , WR-0 , DWRL-0 , BMRL=0 

12  FIXED  J.  - ! 

WL=0 , WR-0 , DWL-0 , DWR-0 

6250  REM  13  E.S.S  J.  - K"  — 

WLR-0 , DWRL=0 , BMRL-0 , SHRL-vk *WR=0 

6260  REM  (in  case  of  a  HOLE) 

6270  REM  11  S.S  e.in  A - 

WR-0 , BMR-0 

12  FIXED  e.in  i - 

WR-0 , DWR-0 

14  F. FIXED  E.in  !  ! - 

WR-0 , SHR-0 

6280  REM  13  E.S.S  e.in  K-- - 

BMR-0 , SHR-vk ( I ) *WR=0 

1 0  FREE  e.in  - 

BMR-0, SHR-0 

6290  REM  Sub.  which  checks  if  the  solution  satisfies  the 
differential  equation 


6300  '  VAR;  n  ,  i  ,  k  ,  sofd  <  i  >,,  so-fed  ( i  ),  f  ( i  ,  k )  df  ( i  ,  k )  , 
d2f  <i  ,k)  ,d3f  (i  ,k>  ,  d4f(i  ,k) 

6310  DE4=RI*RI*RI*RI*D4F < I ,K> ; DE3=2*RI*RI*RI*D3F ( I ,K> ; DE2=- 
RI*RI*  < 1+2*N  C*NC+S0CFD  < I > -RI*RI*S0FD ( I ) ) *D2F ( I , K ) 

6320  DE 1 =R I  *  < 1 +2*NC*NC+S0CFD  ( I >  +R I *R I *S0FD (I ) ) *DF ( I  ,  K  >  - 

(NC*NC* < 4-NC*NC+SGCFD  < I >  +RI*RI*S0FD ( I ) ) >*F(I,K>  : 

DDDD=DE 1 +DE2+DE3+DE4 : IF  DDDD>.0001  THEN  GOTO  6340 
6330  PRINT  .-PRINT  "i  =  " ;  1 5  TAB  (20)  ;  "k  =  "  ;  K;  TAB  (40)  ; 

"de(f  (  I;  "  ,  ";K;  " ) =" ; DE 1 +DE2+DE  3+-DE4  ;  PRINT:  PRINT 
"de (f " ; I; " ,";K; DE1+DE2+DE3+DE4: RETURN 
6340  PRINT; PRINT  "  i  =  "  5  I ; TAB (20) ; " k  =  " ; K; TAB (40) ; 

"da(f ("; 1 5 " , ";K; " ) =" ; DE1+DE2+DE  3+DE4 
:  PRINT:  BEEP:  BEEP:  PRINT  •'<<<<<<  you  have  a  mistake 
>>»>>>»"  :  BEEP:  RETURN:  STOP 

6350  REM  printout  of  then  results  (PRIMARY  and  BUCKLING) 

6360  REM  OPEN  "b : out . dat "  FOR  OUTPUT  AS  #1 
6370  REM  PE (0) =1 : E (2) =1 : NPL=2: RES (0) =0: RES ( 1 ) =2: 

I  BOW  ( 0 )  =  1  :  I BCW  ( 1 )  =  1 0 :  I  BCWG$  (0)=M  •••,':IBCWG*(1>=" 

"  :  PE ( 1 ) = 1 

6380  BEEP: BEEP: BEEP: BEEP 

6390  PR  I  NT  SPRINT: PR  I NT "  <  < < < < <  < < <  RESULT 

> >  > > > > > > > > > > > > > " : PR I NT* PRINT 

6400  PRINT  # 1 , : PR I NT  #1 PRINT  41, "  < < < < < < < < <  RESULT 

> > > > > > > > > > > >»>":  PRINT  41,:  PRINT  41, 

6410  PRINT : PRINT"  in-plane  b.c.  " ; 

6420  PRINT  #1,: PRINT  41,"  in-plane  b.c.  " ; 

6430  SI*-" - " 1 S2-$="=====" : SM$="  " : SF*="  !":SE$=" 

" : SUM$=" ":3UMR*=" " 

6440  SDE* - "  " SUMLF$= 1 *■"->" : SF2*= " < - " 1 SUMRF*= " " 

6450  IF  E(NPL)-0  THEN  NPLMAX=NPL-1  ELSE  NFLMAX=NPL 
6460  FOR  1  =  1  TO  NF'LMAX 

6470  IF  IRES  < 1-1 ) =0  THEN  SUM$=SUM* +SE*  ELSE  IF  IRES(I-1)=1 
THEN  SUM*«SUM*+SF*  ELSE  SUM$=5UM*+SM$ 

6480  IF  PE ( I - 1 ) =0  THEN  SUMLF£=SUMLF$+SE£  ELSE  IF  PE(I-1)>0 
THEN  5UMLF*=SIJMLF*+SF1*  ELSE  SUMLr  ;£=SUMLF$+SF2$ 

6490  SUMLF*=SUMLF*+3DE$ 

6500  IF  1/2- I NT ( 1/2) =0  THEN  SUM*=SUM$+S2*  ELSE  SUM$=SUM$+S1$ 
6510  NEXT  I 

6520  IF  E ( NPL ) *0  THEN  GOTO  6530  ELSE  GOTO  6550 
6530  I=NPL: IF  IRES(I-1)=0  THEN  SUM*»SUM*+SE$+SDE$  ELSE  IF 
IRES  < 1-1 ) =1  THEN  SUM*»SUM*+SF*+SDE*  ELSE 
SUM*=SUM$+SM$+SDE$ 

6540  IF  PE ( 1-1 ) =0  THEN  SUMLFS=SUMLF*+SE$+SDE*  ELSE  IF  PE(I- 
1)>0  THEN  SUMLF$=SUMLF*+SF1*+SDE$  ELSE 
SUMLF*=SUMLF*+SF2$+SDE$ 

6550  IL=LEN  (SUMO  .-FOR  I  =  IL  TO  1  STEP  -1: 

SUMR*=SUMR*+MID$(SUM$ , I  ,1) 

6560  3FF$=MID*(SUMLF$, I , 1 ) : IF  SFF$=">"  THEN  5UMRF$=SUMRF*+"< " 
ELSE  IF  SFF*="<"  THEN  SUMRF*>SUMRF >"  ELSE 
SUMRF*=SUMRF4+SFF* 

6570  NEXT  I 

6580  PRINT  TAB (20) ;SUMOSUMRO PRINT  :PRINT  "  Forces: 

TAB (20) ; SUMLFOSUMRF* 


6590  PRINT  #1 , TAB (20) ; SUM*+SUMR*s PRINT  #1,: PR I NT  #1,"  Forces: 

" ;  TAB  <20 ) ; SUMLF$+SUMRF* 

6600  PRINT s PRINT"  Vertical  B.C  " ; 

6610  PRINT  #1,  -.PRINT  #1,"  Vertical  B.C  “ ; 

6620  Sl*=" - "  i  S2#=  "  =====  "  :  SUM**'*  "  s  SUMR*="  " 

6630  SDE#="  . . SUMLM*=" 

" : SUMRJ*=" " s  SUMRM*=" " : SDE2$="  " 

6640  IF  E ( NPL ) =0  THEN  NPLMAX=NPL-1  ELSE  NPLMAX=NPL 
6650  FOR  1=1  TO  NPLMAX 
6660  SUM$=SUM*+IBCWG*< 1-1 > 

6670  SUMLJ#=SUMLJ*+STR*  < I  —  1 ) +SDE$s  SUMLM$=SUMLM*+STR* < I >  s I F 
I < NPLMAX  THEN  SUMLM*=SUNLM*+SDE$  ELSE 
SUMLM*=SUMLM$+SDE2* 

6680  IF  I /2-INT  < I /2>  =0  THEN  SUM#-SUM#+S2#  ELSE  SUM#-SUM#+S1# 
6690  NEXT  I 

6700  IF  E (NPL) =0  THEN  GOTO  6710  ELSE  GOTO  6720 
6710  I =NPL : SUM#«SUM#+ I BCWGf < I - 1 ) + " 

" : SUMLJ#=SUMLJ$-t-STR*  < 1-1 ) +SDE#i SUMLM*=SUMLM$+" 

6720  I L=LEN  <  SUM# )  :  FOR  I  =  IL  TO  1  STEP  - 
1 : SUMR#=SUMR#+M I D* <  SUM# ,1,1) 

6730  SUMRJ#=SUMRJ#+MID#  (SUMLJ#,  I  ,  1)  .-NEXT  I 
6740  IL»LEN(SUMLM#) : FOR  I=IL  TO  1  STEP  - 

1 s  SUMRM#=SUMRM#+M I D# ( SUMLM# , I , 1 ) : NEX  T  I 
6750  PRINT  T AB  < 20 ) ; SUM# +SUMR# : PR  I NT 

TAB  <20) sSUMLJ#+SUMRJ#sPRINT  TAB (20) 5 SUMLM4+SUMRM* 

6760  PRINT  #1 , TAB (20) s SUM#+SUMR#s PRINT 

#1 , TAB ( 20 ) s SUMLJ #+SUMRJ# : PR I NT  #1 , TAB (20) 5 SUMLM*+SUMRM* 
6770  IF  IPRI#="no"  THEN  GOTO  6780  ELSE  IPRI#="no"  1  RETURN 
6780  PR I NT : PR I NT : PR I NT "  <<<<<<<<<  data  and  results 
> > >  >  >  >  > >  > >  > > "sPRINT 

6790  PRINT  #1, SPRINT  #1, sPRINT  #1,"  <<<<<<<<<  data  and 
results  > > > > > > > > > > > > ".'PRINT  #  1  , 

6800  PR I NT  s  PR I NT  TAB (IS)?"  <  < <  <  <  “ 5  NAMEC# ; " 

>>>>»"  1  PRINT 

6810  PR I NT # 1 , : PR I NT  #l,TAB<15>s"  <<<<<  " ; NAMEC$ ; " 

>>»>>"  1  PRINT  #1, 

6820  SS3#= " ### . ### " s  SS2#= " ## . ## ” 

6830  PRINT  TAB (30) 5 "nc (no.  of  mode) ="5  NC 

6840  PRINT  TAB (30) j "po (poison  rat .)  =  "; 1  PRINT  USING  882# 1  POO 
6850  PRINT  TAB (30) 5 "th (thickness  >="; 1  PRINT  USING  SS2#sT(l> 
6860  PRINT  #1 ,TAB (30) 5 "nc (no.  of  mode)=";NC 
6870  PRINT  #1 ,TAB(30) 5 "po(poison  rat . ) =" 5 : PRINT  #1, USING 
SS2#j  POO 

6880  PRINT  #1 , TAB (30) 5 "th (thickness  )=" 5 sPRINT  #1, USING 
SS2#s  T ( 1 ) 

6890  PRINT  sPRINT  sPRINT  TAB (20) 5"  Buckling  Force  (Pcr*D/R/R> 

"SPRINT  TAB (19) ; " - 

"sPRINT 

6900  PRINT  #l,s PRINT  #1, sPRINT  # 1 , TAB ( 20) 5 "  Buckling  Force 

(Pcr*D/R/R>  "sPRINT  41 , TAB ( IS) • " - 

- "-PRINT  #1i 

6910  IF  IYEINTE=1  THEN  I  YE INTE=Os RETURN 


6920  PRINT  TAB<1> 3 "rl/rO  \e2/el " 5  TAB ( 15) ; : PRINT  USING 
SS3$; E22 (1 )  5  s PRINT  TAB < 25) ; : PRINT  USING 
SS3*sE22(2> PRINT  TAB (35) 3 s PRINT  USING 
SS3*;E22(3) ;  SPRINT  TAB (45) s: PRINT  USING  SS3$;E22(4>s 
6930  PRINT  TAB (55) \ s  PRINT  USING  SS3*s E22 (5) 5 1  PRINT 

TAB (65) PRINT  USING  SS3$; E22 (6) ; s PRINT  TAB (75) PRINT 
USING  SS3$sE22(7> SPRINT 

6940  FOR  I R 1 0= 1  TO  9s  RR I® I RIO/ 10: PRINT  TAB (4) PRINT  USING 
SS2*;RRI; 

6950  FOR  IE21=1  TO  7:PRINT  TAB  (3+IE21*  10)  PRINT  USING 
SS3* ;  PCRR  ( I E2 1 ,  IR10)*R<0)  *R(0>  ; 

6960  NEXT  IE21 : NEXT  IR10 

6970  PRINT  #1 ,  TABU)  ;  "rl/rO  \e2/el TAB  ( 15)  ::  PRINT  #1,  USING 
SS3$ : E22 ( 1 ) ; : PR  I NT  #1 , TAB (25) ; : PRINT  #1, USING 
SS3$;  E22  (2) :PRINT  #1  , TAB  (35)  PRINT  #1,  USING 
SS3S; E22 (3) j SPRINT  #1 , TAB (45) ; s PRINT  #1, USING 
SS3#s  E22 (4) : 

6980  PRINT  #1 , TAB (55) s SPRINT  #1, USING  SS3*| E22 (5) s s PRINT 
#1 ,TAB(65) 5  SPRINT  #1, USING  SS3$; E22 (6 > 3  s PRINT 
# 1 , T AB ( 75 ) ( : PR I NT  #1, USING  SS3* 5 E22 (7) : PRINT 
6990  FOR  I R 10=0  TO  9s RRI=IR10/ 10: PRINT  #1 , TAB (4) 5 : PRINT 
#1, USING  SS2$:RRI; 

7000  FOR  IE21  =  1  TO  7s  PRINT  #1 , TAB <3+IE21*10) s s PRINT  #1 , USING 
SS3$; PCRR ( IE21 , IR10) *R (0) *R (0) ; 

7010  NEXT  IE21 s  NEXT  IR10 
7020  RETURN 

7030  REM  Sub.  far  printout  of  normal  mode 
7040  SS1*="##. ##" s  SS2*="###. ###" s  SS3*="####. ###" 

7050  PRINT  SPRINT  TAB (25);"  Normal  Mode  #"3NC;" 

7060  PRINT  TAB  (25);" - " 

7070  PR I NT  SPRINT SPRINT 

TAB  (28)  :  "Pcr  =  "  ;  PCRR  ( IR10,  IE21 > *R (0) *R (0) / 1 ;  "  DO /r  (0)  "-2" 
7080  PCRR 1 =PCRR ( I R 1 0 , I E2 1 ) sPRINTsPRINT 
7090  PRINT  TAB (3) ; "rr " ; TAB ( 14) ; "Nrr " ; TAB (24 ) ; "Noo" ; 

TAB (33) ; "U (Hor . > " ; TAB  (43); 

" W ( Ver . ) " ; TAB (53) ; "Dw/dr " ; TAB (63) ; "B.mom"; 

TAB ( 73 ) ; " Shear " s  PR  I  NT 
7100  FOR  1  =  1  TO  NPLI s  FOR  J=0  TO  10 

7110  PRINT  TAB  < 1 ) ; : PRINT  USING  SS1$; RR ( I , J ) /R (0) ; s PRINT 
TABUO)  ;  SPRINT  USING  SS2f  l  PCRR1  *SRR  ( I  ,  J  )  ;  :  PRINT 
TAB (20) ; SPRINT  USING  SS2*; SOO ( I , J ) *PCRR 1 ; : PRINT 
TAB (30) : SPRINT  USING  SS2$; UU ( I , J) *PCRR1 ; 

7120  PRINT  TAB (40) ; SPRINT  USING  S33#; WW < I , J ) /WMAX ; s PRINT 
TAB<50) 1  SPRINT  USING  SS3$s  DWW ( I , J) /WMAXs  sPRINT 
TAB (60) s SPRINT  USING  883* S BMW ( I , J ) /WMAX s s 
PRINT  TAB (70) s sPRINT  USING  SS3$s SHW ( I , J > /WMAX : NEXT 
J  s  NEXT  I 

7130  IF  E  ( NPL )  =0  THEN  GOTO  7160  ELSE  I=NF‘Ls  FOR  J=0  TO  5 
7140  PRINT  TABU)  s  sPRINT  USING  SSl*s  RR  ( I  ,  J  >  /R  (0)  5  :  PRINT 
TABUO)  s  sPRINT  USING  SS2*s  PCRR1  *SRR  ( I  ,  J  )  3  s  PRINT 
TAB  (20)  s  sPRINT  USING  SS2*;  SOO  ( I  ,  J  )  *F'CRR1  s  :  PRINT 
TAB (30) 5 sPRINT  USING  SS2$; UU ( I , J ) *PCRR1 ; 


7150  PRINT  TAB (40) ; SPRINT  USING  SS3*j WW ( I , J > /WMAX 5 s PRINT 
TAB (50) ; SPRINT  USING  SS3#; DWW < I , J ) /WMAX ; : PRINT 
TAB (60) ; a  PRINT  USING  SS3$; BMW < I , J ) /WMAX; s 
PRINT  TAB(70) ; SPRINT  USING  SS3* ; SHW < I , J > /WMAX s NEXT  J 
7160  PRINT  #1, sPRINT  #1 , TAB (25) ; "  Normal  Mode  #"5NC5" 

7170  PRINT  #1  , TAB  (25)  ;  " - " 

7180  PRINT  #1, sPRINT  #1, sPRINT 

#1 , TAB (28) 5 "Pcr=" ; PCRR ( I  RIO , IE2 1 ) *R ( 0) *R  < 0) / 1 5 " 

DO/r  (0)  '"•2" 

7190  PRINT  #1, sPRINT  #1, 

7200  PRINT  #1 , TAB (3)  ;  "rr "  ;  TAB ( 14) ; "Nrr " ; TAB (24) ; "Noo"; 

TAB (33) ; "U(Hor. ) "s  TAB (43) ; "W(Ver. > “ ; TAB (53) ; 

"Dw/dr " : TAB ( 63) s "B.  mom" 5  TAB (73) s "Shear " s PRINT  #1 , 

7210  FOR  1*1  TO  NPLI s  FOR  J=0  TO  10 
7220  PRINT  #1 ,TAB(1> ; sPRINT  #1, USING 

SSlf; RR ( I , J ) /R (0) ; 5  PRINT  #1 , TAB ( 10) ; s PRINT  #1, USING 
SS2$$PCRR1*SRR(I ,J) 5 sPRINT  #1 , TAB (20) ; s PRINT  #1, USING 
SS2$ ; SOO ( I , J ) *PCRR 1 ; s  PR I NT  #1 , TAB (30) ; s PRINT  #1, USING 
SS2$ ; UU  < I , J  >  *PCRR 1 ; 

7230  PRINT  #1 ,TAB(40) 5 sPRINT  #1, USING 

SS3$sWW ( I ,J) /WMAX; : PRINT  #1 , TAB (50) 5 8  PR  I NT  #1, USING 
SS3$j DWW ( I , J) /WMAXs s PRINT  #1 , TAB (60 ) 5 s PR INT  #1 , USING 
SS3$a  BMW ( I , J ) /WMAX ;  s PRINT  #1 , TAB (70) 5 s PRINT  #1, USING 
SS3*j SHW ( I, J) /WMAXs NEXT  JsNEXT  I 
7240  IF  E (NPL) =0  THEN  GOTO  7270  ELSE  I=NPLsFOR  J=0  TO  5 
7250  PRINT  #1 , TAB (1) 5 sPRINT  #1, USING 

SS1$  j  RR  < I , J > /R (0) 5  s PRINT  #1 , TAB ( 10) s s PRINT  #1, USING 
SS2#?PCRR1*SRR< I ,J> ; sPRINT  #1 , TAB (20) ; s PRINT  #1, USING 
SS2$sS00< I ,J)*PCRRls sPRINT  #1 , TAB (30) ; s PRINT  #1, USING 
SS2* ; UU ( I , J ) *PCRR 1 ; 

7260  PRINT  #1 , TAB (40) 5 sPRINT  #1, USING 

SS3*;WW(I ,J) /WMAXs sPRINT  #1 , TAB (50) ; s PRINT  #1, USING 
SS3$; DWW( I ,J) /WMAX; sPRINT  #1 , TAB (60) 5 s PRINT  #1, USING 
SS3*; BMW ( I , J ) /WMAX ; s  PRINT  #1 , TAB (70) 5 s PRINT  #1, USING 
SS3$ ; SHW (1,0) /WMAX  s  NEXT  J 
7270  RETURN 


